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ABSTRACT 

We model the temperature and chemical structure of molecular clouds as a 
function of depth into the cloud, assuming a cloud of constant density n illu- 
minated by an external FUV (6 eV < hv < 13.6 eV) flux Go (scaling factor in 
multiples of the local interstellar field). Extending previous photodissociation 
region models, we include the freezing of species, simple grain surface chem- 
istry, and desorption (including FUV photodesorption) of ices. We also treat the 
opaque cloud interior with time- dependent chemistry. Here, under certain con- 
ditions, gas phase elemental oxygen freezes out as water ice and the elemental 
C/O abundance ratio can exceed unity, leading to complex carbon chemistry. 
Gas phase H 2 and O2 peak in abundance at intermediate depth into the cloud, 
roughly Ay ~ 3 — 8 from the surface, the depth proportional to ln(Go/n). Closer 
to the surface, molecules are photodissociated. Deeper into the cloud, molecules 
freeze to grain surfaces. At intermediate depths photodissociation rates are at- 
tenuated by dust extinction, but photodesorption prevents total freezeout. For 
Go < 500, abundances of H 2 and 2 peak at values ~ 10~ 7 , producing columns 
~ 10 15 cm -2 , independent of G and n. The peak abundances depend primarily 
on the product of the photodesorption yield of water ice and the grain surface 
area per H nucleus. At higher values of Go, thermal desorption of O atoms from 
grains enhances the gas phase H 2 peak abundance and column slightly, whereas 
the gas phase 2 peak abundance rises to ~ 10~ 5 and the column to ~ 2 x 10 16 
cm -2 . We present simple analytic equations for the abundances as a function of 
depth which clarify the dependence on parameters. The models are applied to 
observations of H 2 0, 2 , and water ice in a number of sources, including B68, 
NGC 2024, and p Oph. 
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1. INTRODUCTION 

Oxygen is the third most abundant element in the universe, after hydrogen and helium, 
so that a basic knowledge of oxygen chemistry in molecular clouds is essential in order to 
understand the chemical structure, thermal balance, and diagnostic line emission from star- 
forming molecular gas in galaxies. Early gas-phase chemical models (e.g., Langer & Graedel 
1989, Millar 1990, Bergin et al 1998) predicted large abundances of H 2 (~ 1(T 6 ) and 2 
(~ 1CT 5 — 10~ 4 ) relative to hydrogen nuclei in molecular gas well shielded from far ultraviolet 
(FUV, 6 eV < hu < 13.6 eV) photons. If gas phase H 2 and 2 were that abundant, they 
would be important coolants for dense gas (Goldsmith & Langer 1978, Hollenbach 1988, 
Neufeld, Lepp, & Melnick 1995). However, the Submillimeter Wave Astronomy Satellite 
(SWAS) made pointed observations in low energy transitions of ortho-H 2 (the lio — loi 
transition at 557 GHz) and 2 (the 33 — 1 2 transition at 487 GHz) toward numerous dense 
(but unshocked) molecular cores and determined that the line-of-sight-averaged and beam- 
averaged (SWAS beam ~ 4') abundance of H 2 is of order < 3 x 1(T 8 (Snell et al 2000) and 
that of 2 is < 10~ 6 (Goldsmith et al 2000). More recent observations by the Odin mission 
set more stringent upper limits on 2 , < 10~ 7 (Pagani et al 2003), with a reported detection 
at the ~ 2.5 x 10~ 8 level in p Oph (Larsson et al 2007). Although the water abundance 
derived from the observed water emission depends inversely on the gas density, and therefore 
is somewhat uncertain, understanding the two order of magnitude discrepancy between the 
gas phase chemical models and the observations is essential to astrochemistry and to the 
basic understanding of the physics of molecular clouds. 

Previous attempts to explain the low abundance of H 2 and 2 observed by SWAS 
showed that time dependent gas-phase chemistry by itself would not be sufficient (Bergin 
et al 1998, 2000). Starting from atomic gas, a dense (n(H 2 ) ~ 10 5 cm~ 3 ) cloud only took 
> 10 4 years to reach H 2 abundances ~ 10~ 6 , close to the final steady-state values and much 
greater than observed. 

The best previous explanation involved time-dependent chemistry linked with the freeze- 
out of oxygen species on grain surfaces and the formation of substantial water ice mantles 
on grains (Bergin et al 2000). In these models, all the oxygen in the molecular gas which is 
not tied up in CO adsorbs quickly (in ~ 10 5 years for a gas density of ~ 10 4 cm" 3 ) to grain 
surfaces, forms water ice and remains stuck on the grain as an ice mantle. As a result, the 
gas phase H 2 abundance drops from ~ 10~ 6 to ~ 10~ 8 . The grains are assumed too warm 
for CO to freeze as a CO ice mantle, so that for a period from about 10 5 years to about 
3 x 10 6 years, the gas phase H 2 abundance remains at the ~ 10" 8 level, fed by the slow 
dissociation of CO into O and the subsequent reaction of some of this O with H^, which 
then ultimately forms gas-phase H 2 0. The slow dissociation of CO is driven by cosmic ray 
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ionization of He; the resultant He + reacts with CO to form O and C + . After ~ 3 x 10 6 years, 
even the gas phase CO abundance drops significantly as the dissociation process bleeds away 
the oxygen which ends up as water ice on grain surfaces. Therefore, after about 3 x 10 6 years, 
the gas phase H 2 also drops from ~ 10~ 8 to < 10~ 9 . This previous scenario explained the 
observed H2O emission as arising from the central, opaque regions of the cloud, where the 
abundance has dropped to the observed values but have not had time to reach the extremely 
low steady state values. The model relied, then, on a "tuning" of molecular cloud timescales 
so that they are long enough for the freeze-out of existing gas phase oxygen not in CO onto 
grains, but not so long that the CO is broken down and the resultant O converted to water 
ice, which would cause the gas phase H 2 abundances to drop below the observed values^] 
The model also relied on CO not freezing out on grains in the opaque cloud centers. 

To add to the problems of modeling the SWAS data, the SWAS team observed a number 
of strong submillimeter continuum sources such as SgrB, W49 and W51 and found the 
557 GHz line of H 2 in absorption, as the continuum passed through translucent clouds 
(A v ~ 1 - 5) along the line of sight (Neufeld et al 2002, Neufeld et al 2003, Plume et al 
2004). The absorption measurement provided an even better estimate of the H2O column, 
ATi 2 0, through these clouds because the absorption line strengths are only proportional to 
N(H 2 0), whereas the emission line strengths are proportional to n(H 2 ) N (H 2 0)e~ 27K ^ T since 
the line is subthermal, effectively optically thin, and lies AE/k = 27 K above the ground 
state. Therefore, to obtain columns from emission line observations requires the separate 
knowledge of both the gas density and the gas temperature. The absorption measurements 
showed column-averaged abundances of H 2 of ~ 10~ 7 — 10~ 6 with respect to H 2 , using 
observed CO columns and multiplying by 10 4 to obtain H 2 columns. In the context of the 
time-dependent models with freezeout, the difference between the abundances measured in 
absorption compared with those measured in emission was attributed to the presumed lower 
gas densities in the absorbing clouds compared to the emitting clouds. Because the freezeout 
time depends on n _1 , it was assumed that the lower density clouds had not had time to freeze 
as much oxygen from the gas phase, so that gas phase H 2 abundances were higher. 

A variation of this model is that of Spaans & van Dishoeck (2001) and Bergin et al 
(2003), where it was noticed that the water emission seemed to trace the photodissociation 
region (PDR, see Hollenbach & Tielens 1999), which lies on the surface (Ay<5) of the 
molecular cloud. A two component model was invoked in which the water froze out as ice 
in the dense clumps, but remained relatively undepleted in the low density interclump gas 



'■For another similar model which relies on time-dependent effects, coupled with grain freezeout and grain 
surface processes, to explain the low O2 and H2O abundances observed by SWAS, see Roberts & Herbst 
(2002). 
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because of the longer timescales for freezeout. Thus, the average gas phase water abundance 
was derived from a mix of heavily depleted and undepleted gas. 

We propose in this paper a new model of the H 2 and 2 chemistry in a cloud. In this 
model, we assume that the molecular cloud lifetime is sufficiently long to allow freezeout^ 
and reduce the gas phase H 2 and O2 abundances to very low values in the central, opaque 
regions of the cloud. The key to understanding the H2O emission and the O2 upper limits 
observed by SWAS and Odin is to model the spatially- dependent H 2 and 2 abundances 
through each cloud. At the surface of the cloud the gas-phase H 2 and O2 abundances are 
very low because of the photodissociation by the ISRF or by the FUV flux from nearby OB 
stars. Near the cloud surface the dust grains have little water ice because of photodesorption 
by the FUV field. Deeper into the cloud, the attenuation of the FUV field leads to a rapid 
rise of the gas-phase H 2 and O2 abundances, which peak at a hydrogen nucleus column 
Nf ~ 10 2L5 — 10 22 cm -2 (or depths Ayj ~ several into the cloud) and plateau at this value 
for A Ay ~ several beyond Ayf, insensitive to the gas density n and the incident FUV flux 
Go (scaling factor in multiples of the avearage local interstellar radiation field). At these 
intermediate depths, photodesorption of some of the water ice by FUV photons keeps the gas 
phase water abundance high (the "/" in Nf and Ayf signifies the onset of water ice freezeout, 
as will be discussed in §2.5 and §3.3). The FUV is strong enough to keep some of the ice off 
the grains, but, due to efficient water ice formation followed by photodesorption, it is not 
strong enough to dissociate all gas-phase water. At greater depths, gas phase reactions other 
than photodissociation begin to dominate gas phase H 2 destruction, and the steady state 
gas phase H 2 and O2 abundances plummet as all the gas phase elemental O is converted 
to water ice. Thus, the overall structure of a molecular cloud has three subregions: at the 
surface is a highly "photodissociated region", at intermediate depths is the "photodesorbed 
region" , and deep in the cloud is the "freezeout region" . 

In this new model, the H 2 and O2 emission mainly arise from the layer of high abun- 
dance gas in the plateau which starts at Nf and extends somewhat beyond, so that the 
emission is a (deep) "surface" process rather than a "volume" process. For clouds with 
columns N > Nf, the H2O and O2 emission becomes independent of cloud column, as- 
suming the incident FUV field and gas density are fixed. The model gives column- averaged 
abundances through the dense cores of < 10~ 8 for H 2 and 2 for low values of G < 500, 
but the local abundances peak at values which are at least an order of magnitude larger than 
these values. For cloud columns N > Nf, the average H 2 abundance scales as N' 1 . 



2 Indeed, the freezeout of volatile molecules such as CO is now commonly observed in the dense, cold 
regions of molecular clouds (e.g., Bergin & Tafalla 2007) 



- 5 - 



The implications of this scenario for molecular clouds may be much broader than this 
particular model of H 2 and O2 chemistry in clouds. Other molecules such as CO, CS, CN, 
and HCN require a spatial model of their distribution in a cloud, with photodissociation 
and photodesorption the dominant processes near the cloud surface, and the freezing out 
of the molecules the dominant process deeper into the cloud. In addition, the adsorption 
process and creation of abundant ice mantles changes the relative gas phase abundances of 
the elements. In the case considered in this paper, the C/O ratio in the gas may go from 
0.5 at the cloud surface to unity or greater in the H 2 freezeout region. Such changes in the 
gas phase C/O ratio have major implications on all gas phase chemistry (e.g., Langer et al 
1984, Bergin et al 1997). 

The interesting and important caveat to this relatively simple steady state model is the 
effect of time-dependence in raising abundances of, for example, gas phase H 2 0, 2 , and CO 
above steady state values in the opaque (Ay > 5 — 10) interiors of clouds. One such effect 
is that at low densities, n< 10 3 cm~ 3 , species do not have time to freeze out within cloud 
lifetimes and therefore have much higher gas phase abundances. We have also uncovered a 
new time- dependent process that may elevate the H2O and 2 abundances for times t < 10 7 
years in the freezeout region at very high Ay even at high cloud densities n ~ 10 5 cm~ 3 . 
If the grains are sufficiently cold to freeze out CO (T gr < 20 K or G < 500) at high A v , a 
CO/H 2 ice mix rapidly (t< 10 5-6 years) forms. The steady state solution has very little 
CO ice with most of the O in H 2 ice. However, the time to convert CO ice to H 2 ice 
is very long, and the bottleneck is the cosmic ray desorption of CO from the CO/H 2 ice 
mixture. This desorption provides gas phase CO which then acts as a reservoir to produce 
gas phase H 2 and 2 , until all the oxygen eventually freezes out as water ice (i.e., the same 
mechanism as described in Bergin et al (2000) except the gas phase CO in this case comes 
from the CO ice). Depending on the assumptions regarding the CO desorption process, this 
cosmic ray CO desorption timescale may range from 10 5 to >> 10 7 years. If the timescale 
is roughly 0.1 to 1 times the cloud age, a maximum gas phase H 2 and 2 abundance is 
produced at that time. Although not as large as the peak abundances produced at Ayf, 
these abundances can be significant and can contribute to the total column of these species 
if the cloud has a high total column (but only if the CO desorption time is 0.1 to 1 times 
the cloud age). 

This paper is organized as follows. In §2 we describe the new chemical/thermal model 
of an opaque molecular cloud illuminated by FUV radiation. The major changes imple- 
mented in our older photodissociation region (PDR) models (Kaufman et al 1999) include 
the adsorption of gas species onto grain surfaces, chemical reactions on grain surfaces, and 
the desorption of molecules and atoms from grain surfaces. In §3 we show the results of the 
numerical code as functions of the cloud gas density n, the incident FUV flux Go, and the 
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grain properties. We present a simple analytical model that explains the numerical results 
and how the abundances scale with depth and with the other model parameters. We discuss 
time dependent models for the opaque cloud center. In §4 we apply the numerical model to 
both diffuse clouds and dense clouds. Finally, we summarize our results and conclusions in 
§5. For the convenience of the reader, Table 3 in the Appendix lists the symbols used in the 
paper. 

2. THE CHEMICAL AND THERMAL MODEL OF AN OPAQUE 

MOLECULAR CLOUD 

2.1. Summary of Prior PDR Model and Modifications 

The numerical code we have developed to model the chemical and thermal structure 
of an opaque cloud externally illuminated by FUV flux is based on our previous models 
of photodissociation regions (PDRs), which are described in Tielens & Hollenbach (1985), 
Hollenbach et al (1991), and Kaufman et al (1999,2006). This ID code modeled a constant 
density slab of gas, illuminated from one side by an FUV flux of 1.6 x 10~ 3 Go erg cm~ 2 
s _1 . The unitless parameter Go is defined above in such a way that Go ~ 1 corresponds to 
the average local interstellar radiation field in the FUV band (Habing 1968). The subscript 
"0" indicates that the flux is the incident FUV flux, as opposed to the attenuated FUV flux 
deeper into the slab. The code calculated the steady state chemical abundances and the gas 
temperature from thermal balance as a function of depth into the cloud. It incorporated 
~ 45 chemical species, ~ 300 chemical reactions, and a large number of heating mechanisms 
and cooling processes. The chemical reactions included FUV photoionization and photodis- 
sociation; cosmic ray ionization; neutral-neutral, ion-neutral, and electronic recombination 
reactions, including reactions with charged dust grains; and the formation of H 2 on grain 
surfaces. No other grain surface chemistry was included in these prior models. The main 
route to gas phase H2O and O2 through these gas-phase reactions is schematically shown in 
Fig. 1. It starts with the ionization of H2 by cosmic rays or X-rays, which eventually leads to 
H 3 + that recombines to form gas phase water. It also recombines to form OH which reacts 
with O to form gas phase 2 . The code contained a relatively limited carbon chemistry; we 
included CH and CH 2 , as well as their associated ions, but did not include CH 3 , CH 4 , C 2 H 
or longer carbon chains, nor any carbon molecules that include S or N. 

The code modeled regions which lie at hydrogen column densities N < 2 x 10 22 cm~ 2 (or 
Ay < 10) from the surface of a cloud. Therefore, it applied not only to the photodissociated 
surface region, where gas phase hydrogen and oxygen are nearly entirely atomic and where 
gas phase carbon is mostly C + , but also to regions deeper into the molecular cloud where 
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hydrogen is in H2 and carbon is in CO molecules. Even in these molecular regions, the 
attenuated FUV field can play a significant role in photodissociating H 2 and O2, and in 
heating the gas. 

The dust opacity is normalized to give the observed visual extinction per unit H nucleus 
column density that is observed in the diffuse ISM (e.g., Savage & Mathis 1979, Mathis 
1990). The extinction in the FUV, assumed to go as e~ bxAy , is taken from PDR models 
described most recently in Kaufman et al (1999, 2006), and mainly derived from the work of 
Roberge et al (1991). Here b\ is the scale factor which expresses the increase in extinction 
as one moves from the visual to shorter wavelengths. However, the FUV extinction law in 
molecular clouds is not well known. We discuss in §3.2 how our results would vary if b\ 
varies. 

In the following subsections (§2.2 - §2.7) we describe the modifications to the basic PDR 
code which we have made to properly model the gas phase H2O and O2 abundances from 
the surfaces of molecular clouds to the deep interiors, where the gas phase elemental oxygen 
freezes out as water ice mantles on grain surfaces. The modifications include the freezeout 
of gas phase species onto grains, the formation of OH and H 2 as well as CH, CH 2 , CH 3 
and CH4 on grain surfaces, and various desorption processes of species from grain surfaces. 
Fig. [2] shows the grain surface formation process for water ice. 

One key modification is the application of time dependent chemistry to the deep, opaque 
interior. Here, certain chemical timescales are comparable to cloud lifetimes, and steady 
state chemical networks such as our PDR code do not apply. The time dependent code also 
increases the number of reactions and species we follow, and this particularly improves the 
treatment of the carbon chemistry in the cloud interior. Further discussion is found in §2.3, 
2.4 and 3.7. 



2.2. Dust Properties 

The basic PDR code implicitly assumed a "MRN" (Mathis, Rumpl & Nordsieck 1977) 
grain size distribution n gr (a) oc a -3,5 , where a is the grain radius, extending from very small 
( < 10 A) polycyclic aromatic hydrocarbons (PAHs) to large (~ 0.25 //m) dust grains. With 
this distribution the largest grains provide the bulk of the grain mass, and the intermediate- 
sized grains (~ 100 A) provide the bulk of the FUV extinction. The smallest particles and 
PAHs provide the bulk of the surface area, and dominate the grain photoelectric heating 
process and the formation of H 2 on grain surfaces. The code implicitly assumes an MRN 
distribution in the sense that the equations for attenuation of incident photons by dust, for 
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grain photoelectric heating, for H 2 formation, and for gas-grain collisions and the charge 
transfer that may occur in such collisions all assume such a distribution. 

For the process of the freezeout of oxygen species onto grain surfaces, we also implicitly 
assume an MRN distribution, in the sense that we adopt a cross sectional grain area a H 
appropriate for an MRN distribution. As will be shown in §3, an, which is one fourth the 
surface area available for freezeout, is an important parameter for determining the abundance 
of gas-phase H 2 in our models. Below in §2.5, we show that on is roughly the total cross 
sectional area of grains (per H) that are larger than about 20 A. Smaller grains experience 
single photon heating events which clear them of ice mantles (Leger, Jura, & Omont 1985). 
In addition, cosmic rays may clear small grains of ices (Herbst & Cuppen 2006). In an 
MRN distribution, the cross sectional area of grains with radii larger than about 20 A is 
ff fl ~2x 10~ 21 cm 2 per H. 

We note, however, that the cross sectional area is somewhat uncertain. It could be 
smaller if grains coagulate inside of dense molecular clouds. It could be larger due to ice 
buildup. We note that once all oxygen-based ices freeze out onto such an MRN distribution, 
it forms a constant mantle thickness of Aa ~ 50 A regardless of grain size. While this is 
insignificant to the large (~ 0.1 jum) grains which contain most of the mass, it represents 
a very large change in radius and area of the smallest (~ 20 A) grains which provide most 
of the grain surface area. The increased surface area of grains at the end of this freezeout 
will proportionally increase subsequent gas-grain interaction rates. Thus, neutralization of 
ions by grains, formation of H 2 on grain surfaces, the freezeout of species, and the cooling 
or heating of the gas by gas-grain collisions may all be enhanced by significant amounts in 
regions where water ice has totally frozen out on grains. However, because of the possibility 
of grain coagulation which reduces grain surface area, and because of the uncertainly in the 
lower size cutoff of effective grains, we have chosen to fix <j H and treat it as a (constrained) 
variable with values ranging from a H = 6 x 10~ 22 to 6 xl0~ 21 cm 2 per H. 



2.3. Adsorption and Time Dependent Chemistry 

Adsorption is the process by which a gas phase atom or molecule hits a grain and sticks 
to the surface due to van der Waals or stronger surface bonding forces. Typically, at the 
low gas temperatures we will consider, T gas < 50 K, the sticking probability is of order unity 
for most species (Burke & Hollenbach 1983). We assume unit sticking coefficients here for 
all species heavier than helium. Thus, the timescale to adsorb (freeze) is the timescale for a 
species to strike a grain surface: 
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t/,i ~ [n gr (JgrVi\ \ (1) 

where Vi is the thermal speed of species i and where the grain number density n gr times 
the grain cross section a gr is the average over the size distribution. As discussed above, 
in an MRN distribution with a lower cutoff size of a ~ 20 A, n gr a gr ~ 2 x 10~ 21 n cm -1 , 
where n is the gas phase hydrogen nucleus density in units of cm -3 . Therefore, tf^ ~ 
8 x 10 4 (m o /m i ) 1 / 2 (10 4 cm- 3 /n)(30K/T) 1 / 2 years, where mo is the mass of atomic O and 
rrii is the mass of the species. Since molecular cloud lifetimes are thought to be 2 x 10 6 to 
2 x 10 7 years, and the molecules reside in regions with n> 10 3 cm -3 , this short timescale 
suggests that molecules should mostly freeze out in molecular clouds unless some process 
desorbs them from grains. Dust vacuums the condensibles quickly! 

Consider the depletion of CO in regions where the grain temperatures are too high to 
directly freeze out CO. CO freezes out when grain temperatures are T 9r <20 K (see §2.4 
below). However, H 2 freezes out when grain temperatures are T gr < 100 K at the densities 
in molecular clouds. Suppose that 20 K< T gr < 100 K, so that gas phase elemental oxygen 
not in CO freezes out as water ice, but CO does not directly freeze. CO is continuously 
destroyed by He + ions created by cosmic rays, producing C + and O. The resultant O has 
two basic paths: it can reform CO, or it can freeze out on grain surfaces as water ice, never to 
return to the gas phase in the absence of desorption mechanisms. The latter process therefore 
slowly depletes the CO gas phase abundance, by lowering the available gas phase elemental 
O. We see this effect in the late time evolution of CO in Bergin et al (2000). One effect, 
not included in Bergin et al, that slows this process is that for grain temperatures > 20 K 
adsorbed atomic O may be rapidly thermally desorbed before reacting with H and therefore 
the production of H 2 ice by grain surface reactions may be substantially curtailed (this will 
be discussed further in §2.5 and §3.2). Nevertheless, gas phase H 2 will form and will freeze 
out, and therefore, given enough time, the elemental O not incorporated in refractory (e.g., 
silicate) grain material is converted to H 2 ice. If elemental C does not evolve to a similarly 
tightly bound carbonaceous ice, there is the potential for the gas phase C/O number ratio to 
exceed unity. This will have consequences for astrochemistry (e.g., Langer et al 1984, Bergin 
et al 1997, Nilsson et al 2000), and we show some exemplary results in §3.5. 

Consider as well the case where T gr < 20 K deep in the cloud (i.e., at Ay > 5; this low 
T gr corresponds to Go <500). In this case, the gas phase CO freezes out with the H 2 0, 
forming a CO/H 2 ice mix. In order for this CO ice to be converted to H 2 ice, the CO 
needs to be desorbed by cosmic rays. It then reacts in the gas phase with He + as described 
above and eventually the freed O goes on to form mostly water ice. The timescale for this 
CO cosmic ray desorption can be quite long, and may dominate the slow evolution of the 
chemistry toward steady state. 
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Unfortunately, this timescale is not well known. First of all, cosmic ray desorption 
rates are not well constrained. But secondly, the timescales depend on the structure of the 
CO/H 2 ice mix, and this structure is not certain. If the molecules are immobile, the CO 
and H 2 ice molecules are initally mixed and each monolayer contains this mix. However, 
there is some evidence that either pure CO ice pockets form, or the CO accretes on top of 
an already formed layer of water ice (Tielens et al 1991, Lee et al 2005, Pontoppidan et al 
2004, 2008). Even if the CO and H 2 are initially mixed and immobile, cosmic rays desorb 
CO much more efficiently than H 2 due to the much smaller binding energy of adsorbed 
CO compared with adsorbed H 2 (Leger, Jura, & Omont 1985). This could lead to surface 
monolayers with very little CO compared to the protected layers underneath the surface. The 
cosmic ray desorption rate is directly proportional to the fraction of surface sites occupied 
by CO. In the examples above, this fraction could range from zero to unity. Thus, the 
desorption timescale could range from large values (>> 10 7 years) to small values (~ 10 5 
years). We compute these timescales in §2.4. 

Because of the moderately long timescale for the depletion of either gas phase or ice 
phase CO via this freezing out of the elemental O in the form of water ice, we apply time 
dependent chemical models (see §3.7) to the cloud interiors which are run for the presumed 
lifetimes (~ 10 7 years) of molecular clouds. Time dependent models are also warranted for 
clouds of relatively low density, n < 10 3 cm -3 , where tf > 2 x 10 6 years (see Eq. 1). 

2.4. Desorption 

To desorb a species from a grain surface requires overcoming the binding energy which 
holds the species to the surface. Table 1 lists the binding energies of the key species we treat. 
Note that water has a very high binding energy and is therefore more resistant to desorption. 

Thermal Desorption. The rate Rtd,i per atom or molecule of thermal desorption from a 
surface can be written: 

R tdtl ~ u i e- Ei ' kT <r, (2) 

where E { is the adsorption binding energy of species i and v { = 1.6 x 10 11 ^ (Ei/k) / {mi/m^) 
s -1 is the vibrational frequency of the species in the surface potential well, k is Boltzmann's 
constant and rrii and vrtu are the mass of species % and hydrogen, respectively. One can find 
the temperature at which a species freezes by equating the flux F t d,i of desorbing molecules 
from the ice surface to the flux of adsorbing molecules from the gas, 
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F t d,i = N S:i R td:i f s>i = 0.25niVi, (3) 

where N 8ji is the number of adsorption sites per cm 2 (jV Sji ~ 10 15 sites cm -2 ), f s>i is the 
fraction of the surface adsorption sites that are occupied by species i, rii is the gas phase 
number density of species % and Vi is its thermal speed. We have assumed sticking probability 
of unity. From these equations, one can calculate the freezing temperature Tf ti of a species 
by solving for T gr . 
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This equation shows that at molecular cloud conditions the freezing temperature Tf^ is about 
0.02Ei/k, or about 100 K for H 2 but only ~ 20 K for CO (see Table 1) We include thermal 
desorption rates from our surfaces as given by these equations. We compute the number of 
layers of an adsorbed species on the grain, and only thermally desorb from the surface layer. 

Photodesorption. The flux F pd i of adsorbed particles leaving a surface due to photodes- 
orption is given by 



F p d,i — YiF FUV f S:i , (6) 

where Yi is the photodesorption yield for species i, averaged over the FUV wavelength band, 
and Fpuy is the incident FUV photon flux on the grain surface. We write the FUV photon 
flux at a depth A v into the cloud as (see Tielens & Hollenbach 1985) 



F FUV = G Q F Q e- 1 -* Av , (7) 

where F ~ 10 8 photons cur 2 s" 1 is approximately the local interstellar flux of 6eV-13.6eV 
photons and G is the commonly used scaling factor (see §2.1). 

Reliable photodesorption yields are often hard to find in the literature, but for the case of 
H 2 a laboratory study has been done for the yield created by Lyman a photons (Westley et 
al 1995a, b). The mass loss at large photon doses is proportional to the photon flux, showing 
that desorption is not due to sublimation caused by absorbed heat. The yield is negligible 
when infrared or optical photons are incident. As an approximation, we will assume the 
average yield in the FUV (912-2000 A) is the same as the Lyman a (1216 A) yield. Most 
of the photodesorption occurs when the incident photon is absorbed in the first two surface 
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monolayers of the water ice (Andersson et al 2006). The value of the yield (~ 1 - 8 x 10~ 3 ) 
is perhaps 10 times smaller than the probability that the incident photon is absorbed in the 
first two surface monolayers, so the photodesorption efficiency from these layers is of order 
10%. We note that this yield from Westley et al is similar to the values that Oberg (private 
communication) has found in preliminary experiments with a broad band FUV source. 

One peculiar result from the Westley et al experiments is that the yield first increases 
with the dose of photons, until it hits a saturation value for that flux. The saturation dose 
is of order 3 xlO 18 photons cm" 2 . In the fields we are considering, Go ~ 1 — 10 5 , this dose is 
achieved in less than 1000 years. Thus, our grains may well be saturated. The interpretation 
of the increase of the photodesorption yield with dose is that photolysed radicals like O, H 
and OH build up on the surface to a saturation value. Chemical reactions of a newly-formed 
radical with a pre-existing radical may lead to the desorption of the product (usually H 2 0) 
and explain the dose dependence of the yield. Adding further weight to this interpretation 
is that the saturated yields decline with decreasing surface temperature (Yh 2 o is 0.008 at 
T gr = 100 K, 0.004 at 75 K, and 0.0035 at 35 K). As temperatures decrease the radicals 
cannot diffuse across the surface as well to react with other radicals, nor can activation 
bonds be as easily overcome. Since presumably radicals like O, H and OH are being created, 
and it is the process of reforming more stable molecules that may help initiate desorption- 
and may kick off neighboring molecules, other molecules may desorb such as OH, H, 2 , H 2 . 
These are indeed observed, but Westley et al report that most of the desorbing molecules 
are water molecules. We note, however, a recent theoretical calculation by Andersson et al 
(2006) which suggests that one might expect OH and H to be major desorption products, 
but confirming that the yield of H 2 photodesorption may still be of order 10~ 3 . We also 
find (see §4.1 below) that adopting a yield for OH+H photodesorption from water ice which 
is twice that of H 2 photodesorption provides better agreement to the OH/H 2 abundance 
ratio measured in translucent clouds. We define Yqh,w as the yield of OH photodesorbing 
from water ice; Yqh, w — 2Yff 2 o- This distinguishes Y Q h,w from Y Q h, which is the yield of 
OH photodesorbing from a surface covered with adsorbed OH molecules. 

The interpretation that radical formation in the surface layer is important in creating 
photodesorption leads one to re-examine the issue of whether grains near the surfaces of 
molecular clouds will indeed be saturated with radicals on their surfaces. While it is true 
that the photon dose is sufficient, this dose was delivered over a long timescale. Westley et 
al (1995b) point out that it is possible in this case that radicals may recombine in the time 
between photon impacts within an area of molecular dimensions (~ 10~ 15 cm 2 ), which might 
result in smaller yields than the saturated ones. 

In light of the uncertainties in the wavelength dependence of the yield over the FUV 
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range, the photodesorption products, and the level of saturation, we adopt Yh 2 o = 10~ 3 and 
Yoh,w = 2 x 10~ 3 as our "standard" values, but we will show models in which these yields are 
varied by a factor of 10 from their standard values, maintaining their ratio of 0.5. We will 
also use the observations of Ay thresholds for water ice formation to constrain the absolute 
values of these yields (§4.2.1), and the observations of OH and H 2 in diffuse clouds to help 
determine their ratio (§4.1). 

Recent laboratory measurements of CO photodesorption imply a yield similar to that 
for H 2 (Oberg et al 2007). Therefore, we adopt Yco — 10~ 3 as our standard value for CO. 
Unless T gr < 20 K, adsorbed CO is anticipated to be rare because it thermally desorbs so 
rapidly. However, in cases with low Go, CO may occupy ~ 1/3 of the surface sites at the 
H 2 freezeout column. 

For the other species there is less or no reliable photodesorption data in the literature. 
Since our grain surface is largely H 2 ice, and the photodesorption is via the production 
of radicals (and possibly the reformation of H 2 from radicals), we assume that the yield 
for O is Yq = 10 -4 (a smaller yield since no reformation possible). The yields for O, OH 
and 2 are less important, because in our models they rarely occupy many of the surface 
sites. The first two are rare because they rapidly react with H to form H 2 0. Table 1 lists 
the photodesorption yields we adopted. 

A complication arises if grains move rapidly from surface regions (higher temperature 
and possible lower density) where H 2 ice can exist to deeper (lower temperature and higher 
density) regions where CO ice forms and then back again to the surface region. The fraction 
of the surface now covered by H 2 ice, f s ,H 2 o, is now nearly zero, since CO covers the 
surface. Therefore, H 2 and OH photodesorption rates are very low. We have ignored this 
effect, assuming a mix of CO and H 2 ice in each adsorbed monolayer^ In this paper the 
region where H 2 and OH photodesorption are important is the intermediate or "plateau" 
region where gas phase H 2 and OH abundances peak. Here, the photodesorption and/or 
thermal desorption timescales of CO are much shorter than the dynamical timescales to 
bring CO-coated ice grains from the deep interior to the plateau region. 

Desorption by Chemical Reaction. There have been some attempts to model the des- 
orption of molecules due to chemical reactions on grain surfaces. D'Hendecourt et al (1982) 
and Tielens & Allamandola (1987) discuss an exploding grain hypothesis in which the UV 
photons produce radicals in the grains, which are relatively immobile until the grain tem- 
perature is increased. At this point, they begin to diffuse and react, and the chemical heat 



3 Except, as discussed later in this section, when we test the sensitivity of our results to CO cosmic ray 
desorption. 
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causes a chain reaction exploding the grain and "desorbing" species in large bursts. Our 
photodesorption is a sort of steady state version of this. O'Neill, Viti, & Williams (2002) 
examine the effects of hydrogenation on grain surfaces, and in some of their models they 
assume that the reaction of a newly adsorbed H atom with radicals on the surface such as 
OH will release enough energy to desorb the resultant hydrogenated molecule. In this case, 
it is the chemical energy carried in by the H atom, and not the photodissociation of the ice 
itself, which leads to desorption. In our standard case we assume that H 2 desorbs on forma- 
tion because of its low adsorption binding energy and its low mass, but that other forming 
molecules do not. However, we ran one case (see §3.5) in which we include this desorption 
process for all forming molecules in our grain surface chemistry. 

Desorption by Cosmic Rays and X-rays. Deep in the cloud, cosmic ray desorption of 
ices becomes important for maintaining trace amounts of O-based and C-based species in 
the gas phase. Cosmic ray desorption rates per surface molecule (i.e., molecule in the top 
monolayer of ice) for various species are listed in Table 1. As discussed above, the desorption 
rates depend on the abundance of a particular adsorbed species in the surface monolayer. 
In our standard steady state models, we assume that if there are a number of species (e.g., 
CO, CH 4 , H 2 0) that are adsorbed, that they are well mixed and each monolayer of ice 
contains the same mix. However, as discussed above, for CO the fraction of the surface layer 
containing CO may range from zero to unity. The timescale for the cosmic ray desorption 
of a species i is given: 



where x(i) is the abundance of the adsorbed ice species and R cr>i is the cosmic ray desorption 
rate per surface atom or molecule and is given in Table 1. The cosmic ray desorption time 
for H 2 Oj ce is long (~ 10 10 years) and is not significant in our models. However, the timescale 
for CO desorption under certain circumstances may be comparable to the age of the cloud, 
and can therefore lead to a source of gas phase CO which can then serve as a source of gas 
phase H 2 and 2 , as discussed above. Substituting from Table 1 for R cr ,co-< assuming an 
initially high abundance of CO ice, x(CO ice ) = 10~ 4 , and taking f s ^co — 0.3, we obtain 
t cr ,co ~ 2.3 x 10 6 years. However, both R cr ,co and f s ,co are uncertain. In the case of 
the fractional surface coverage f s ,co, the cosmic ray desorption process itself, acting on the 
top monolayer, could result in a very low value of f s ,co by essentially purging the surface 
of CO and leaving a fairly pure monolayer of H 2 Oj ce (see §2.3). On the other hand, there 
may be evidence for a non-polar CO ice surface on top of a polar H 2 ice mantle resulting 
m fs,co = 1 (see §2.3). Because of the uncertainly in the product f s ,coRcr,co we will 
examine in §3.7 a range of values of this parameter extending from 10 -15 to 10~ n s _1 . The 
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corresponding timescale to remove all the ice is about 3 x 10 5 to 3 x 10 9 years, ranging from 
less than to much greater than cloud ages. 

Cosmic rays also lead to the production of FUV photons (Gredel et al. 1989) at a level 
equivalent to G ~ 1CT 3 throughout the cloud. We include this flux in determining both the 
photodissociation of gas-phase molecules and the photodesorption of ices. 

Following Leger, Jura, & Omont (1985) we assume that, in general, cosmic ray desorp- 
tion dominates over X-ray desorption and we therefore ignore this process. Of course, in 
regions close to X-ray sources and where the UV is shielded by dust, X-ray desorption will 
become more significant. 



2.5. Timescales Related to Freezeout 

Several timescales are needed to understand the freezing process fully. Assuming that 
the grain size distribution extends to very small grains (or large molecules) whose size is 
of order a few Angstroms, we first compute timescales relevant for single photons to clear 
grains of ice due to transient heating by the photon. The peak temperature that a silicate 
grain of radius a achieves after absorbing a typical ~ 10 eV FUV photon is given (Drapatz 
& Michel 1977, Leger, Jura & Omont 1985) 

/on A\ 13 

T max ~ 100 (-—) K. (9) 
The timescale for this grain to radiate away the energy of this photon is 

100 K\ 5 /20 A x : 



100 ( — I ) seconds. (10) 

max I \ Q> 



On the other hand, the timescale for a single H 2 molecule to thermally evaporate from an 
ice surface on the grain is given 

io 4800 K 

tevap,H 2 ~ 4 X 10 € T ™« seconds, 11 

where AE/k = 4800 K is the binding energy of the H 2 molecule (Fraser et al 2001) 
and 4 x 10~ 13 seconds is the period of one vibration of the molecule in the ice lattice (see 
Eqn 2. and following text). For comparison with the above equations, t evaPi H 2 o — 10 29 
seconds at T max = 50 K, 3 x 10 8 seconds at T max = 100 K, and 30 seconds at T max = 150 
K. Comparison with t r indicates that the ice mantle will evaporate before cooling when 
T max ~ 150 K. Comparison with the equation for T max shows that a > 15 — 20 A is required 
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for grains to be large enough such that single photon events do not clear ice mantles. Thermal 
evaporation of H 2 will also cool the grain. Since the binding energy of an H 2 molecule is 
~ 0.5 eV, ~ 20 H 2 molecules must evaporate to offset the heating due to the absorption 
of a single UV photon. However, we show below that 10 eV photons are absorbed by 20 A 
grains more rapidly than O atoms hit the grains to form ice. Therefore, thermal evaporation 
cannot significantly cool the transiently heated grains compared with radiative cooling, and 
an entire ice mantle never builds up on the grain surface. The conclusion holds that a < 20 
A grains lose their ice in the presence of single photon heating events. In summary, for 
the purposes of modeling oxygen freezeout and ice formation, we assume that the effective 
surface area for ice freezeout is given by the area of grains that are larger than ~ 20 A. 

Several more timescales reveal the fate of an oxygen atom once it sticks to a grain. The 
timescale t gr< o for a given grain of radius a to be hit by a gas phase oxygen atom is given 

10~ 4 \ /10 K\ 1/2 /0.1/mA 2 /10 5 cm"" 



~ 2 ' 5 x 10 {mo)){t^) \tZ-) ) seconds ' (12) 

where x(0) is the gas phase abundance of atomic oxygen relative to hydrogen. Similarly, 
the timescale t gri H for the same grain to be hit by a hydrogen atom from the gas is given 

, , s , /10 K\ 1/2 /0.1/im\ 2 AO 5 cm' 3 \ , . , 

t gr , H ~ 6.0 x lO-^if)- 1 [—J {—^) { n ) seconds ' ( 13 ) 

where x(H) is the abundance of atomic hydrogen. Note that at the column or Ay into 
the cloud where water begins to freeze out on grain surfaces, x(H) is likely less than unity 
because most H is in H 2 . However, if x(H) > 0.25x(O), then H atoms stick more frequently 
than O atoms. Provided photodesorption or thermal desorption does not intervene, and 
assuming that the H atoms are quite mobile on the surface and find the adsorbed O atom, 
O atoms will react with H atoms on the surface to form OH and then H 2 O0 

These timescales need to be compared to the timescale for UV photons to photodesorb 
the O atom and the timescale for O atoms to thermally desorb from grains. The timescale 
for a grain of size a > 100 A to absorb an FUV photon is given by: 



4 In our runs we find only a few cases where x(H) < x(0); this occurs at high Ay for cases with high Go, 
where grains are so warm that O thermally desorbs before forming OH on the grain surface. In these cases, 
there is no buildup of adsorbed O atoms on the surface and therefore no resultant formation of O2 or CO2 
on grain surfaces, as might otherwise be expected (Tielens & Hagen 1982). 
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t 7 ~ 30G 1 e 1 - 8 ^ l^P) seconds. (14) 

The timescale to photodesorb an atom from a surface covered with a covering fraction 
/ S) o of atoms is then, 

i 7 _^ ~ 3 x l^V^V-i (°J£E) 2 _1 seconds, (15) 

where Yo is the photodesorption yield for atomic O. The timescale for thermal desorption 
of an O atom on a grain is given 



t eV ap,o ~ 9 x 10 13 e T max seconds. (16) 

Therefore, the timescale for thermal evaporation of an O atom is 10 22 seconds for a 10 K 
grain, 2 x 10 5 seconds for a 20 K grain, and 0.3 seconds for a 30 K grain. From Eqs (12- 
16) one concludes that if grains are cooler than about 20-25 K and if x(H) > 0.25x(O), an 
adsorbed O atom will react with H to form OH and then H2O before it is thermally desorbed. 
Photodesorption of O never appears to be important. 

In this paper we are mainly interested in the abundance of gas phase H 2 and 2 at 
the point where oxygen species start to freeze out in the form of water ice. This "freezeout 
depth" Ay = Ayf is where the gas phase oxygen abundance plummets and most oxygen is 
incorporated as water ice on grain surfaces. This occurs when the rate at which relatively 
undepleted gas phase oxygen atoms hit and stick to a grain equals the photodesorption rate 
of water molecules from grain surfaces completely covered in ice, or i 7 _ des (Ayf) ~ t gr ,o(^vf)- 
Using Eqs. (12 and 15) and taking Y = Y OHw + Y H , 2Q as the total yield of photodesorbing 
water ice, we find 



A vf ~ 0.56 In 



0.8G 



Y 

1F3 



10 K 



T, 



gas 



1/2 



10 5 cm" 



n 



(17) 



if the term in brackets exceeds unity. If it is less than unity, there is freezeout at the surface 
of the cloud. As can be seen in Eq. (17), freezeout at the surface only occurs for very low 
values of G /n< 10~ 5 cm 3 . 

At Ayf the timescale t 7 / for a grain of size a to absorb an FUV photon can be written, 
substituting Eq. (17) for A V f into Eq. (14), 
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.'10 K\ 1/2 f0.1fim\ 2 /10 5 cm- 3 \ / Y , 
t lf ~ 24 -— — seconds - ( 18 ) 



^gas / V a ) \ n / \ 10 ~ 
As discussed above, t 1 j << t gr> o, so that, at this critical depth, no ice forms on small ( < 20 
A) grains because single photons can thermally evaporate adsorbed O, OH or H 2 before 
thermal radiative emission cools the grain and before another O atom sticks to the grain. 
Although ice may not form on a < 20 A grains, they may still contribute to the formation of 
gas phase OH or H 2 0. Note that once an O atom sticks to a given grain, the timescale t gr ^ 
for an H atom to stick (Eq. 13) may be shorter than the timescales t 1 j for a photon to be 
absorbed. Thus, the O may be transformed into H 2 on the grain surface before a photon 
transiently heats the small grain and thermally evaporates the H 2 to the gas phase. Because 
of these complexities, the possibility of coagulation, and the uncertain nature of small grains, 
we have chosen to fix the cross sectional grain area per H nucleus <jh as constant with depth 
A v in a given model, but to then run models with a range of assumed cr#. 



2.6. Grain Surface Chemistry 

Figure 2 shows the new grain surface oxygen chemistry that we adopt. We maintain the 
old method of treating H 2 formation on grain surfaces (Kaufman et al 1999). The main new 
grain surface chemistry reactions are the surface reactions H + O — >■ OH and H + OH — >• 
H 2 0. Miyauchi et al (2008) discuss and review recent laboratory experiments which suggest 
rapid formation of water ice on grain surfaces. In addition, we include simple grain surface 
carbon chemistry. Similar to Figure 2, we allow the sticking of C to grains to form C ice , 
followed by reactions of surface H atoms to produce CHj ce , CH 2 j ce , CH 3ice , and CH 4 j ce , as 
well as the desorption of all these species. We do not form CO, C0 2 , or methanol on grain 
surfaces in this paper. 

In the model we compare the timescale for O and OH to desorb with the timescale for 
an H atom to hit and stick to the grain surface. If the desorption timescale is shorter, then 
we desorb the radical and the grain surface chemistry does not proceed. However, if the 
desorption time is longer, then we assume that the reaction proceeds. We generally assume 
that the newly formed surface molecule does not desorb from the heat of formation, but must 
wait for another desorption process, or, for a further chemical reaction with H, if possible. 
However, as in O'Neill et al (2002), we also run a case where the heat of formation desorbs 
the newly formed molecule. Implicit in this simple grain chemical code is the assumption 
that the H atom migrates and finds the O or OH before it desorbs or before another H atom 
adsorbs to possibly form H 2 . We note that, in general, it is rare to find even a single O or 
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OH on the grain surface, because desorption or the reaction with H removes them before 
another O or OH sticks to the surface. 



2.7. Ion-Dipole Reactions 

A final change to our chemical code is the inclusion of enhanced rate coefficients for 
the reactions between atomic or molecular ions with molecules of large dipole moment, such 
as OH and H 2 0. We have used the UMIST rates (e.g., Woodall et al 2007) and the rates 
quoted in Adams et al (1985) and in Herbst & Leung (1986). These rates are larger than 
the rates without the dipole interactions by a factor ~ 10 at T = 300 K and, given their 
T -1 / 2 dependence, increase further at low temperatures. 



3. MODEL RESULTS 

In this section, we present computational results of our new model. In order to explore 
the sensitivity of the freeze-out processes to various input parameters, we first define a 
standard case. We then present a number of models in which the gas density, FUV field, 
grain propertie s and photoelectric yie ld are varied. For each of these runs, we use the steady- 



state model of iKaufman et all (119991 ) modified to include grain surface reactions, ion-dipole 
reactions, and desorption mechanisms as detailed in §2. Unlike our previous PDR models 
which were run to a depth Ay = 10, we often run our new models to Ay = 20 to ensure 
that we are well beyond the regions where photodesorption plays a role. At the end of this 
section, we discuss the effects of time-dependence. 



3.1. Standard Case 



Most ortho-H20 detections in ambient (non-shocked) gas are located in Giant Molecular 
Clouds (GMCs), and, in particular, in the massive cores of these GMCs which are illuminated 
by nearby OB stars. Ortho-H 2 has not yet been detected toward quiescent "dark clouds" 
such as Taurus. Therefore, for our standard case, we choose a gas density n = 10 4 cm -3 and 
an FUV fi eld strength Go = 10 2 - We adopt the same gas ph ase abundances as in the PDR 



models of IKaufman et all (119991 ) and IKaufman et al.l (120061 ). Details are given in Tables 1 
and 2. 



In the standard case, the gas (grain) temperature varies from T = 120 K (T gr = 31 K) 
at the surface, to T = 22 K (T gr = 15 K) at the water peak plateau, to T = 16 K (T gr = 15 
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K) for Ay = 10. The abundances of O-bearing species in the standard model are shown in 
Fig. [31 while those of C-bearing species are shown in Fig. HI Here, the basic features of the 
results of freeze-out and desorption are clearly demonstrated. 

The Photodissociation Layer. At the surface, the dominant species are precisely those 
found in gas-phase PDR models. Atomic oxygen is the most abundant O-bearing species, 
while the transition from C + to CO occurs once significant shielding occurs at Ay ~ 1.5. 
Near the surface, all of the oxygen not in CO is atomic, and x(0) > x(CO). At the low dust 
temperature at the surface, thermal desorption of H 2 is insignificant and the beginnings 
of an ice layer (less than a monolayer) form despite rapid photodesorption by the FUV. In 
steady-state, photodesorption of this ice as well as gas phase chemistry maintains a small 
abundance of gas-phase H 2 0. 

The Photodesorbed Layer. By Ay ~ 2.8, enough water has been deposited on grain 
surfaces to form a monolayer, that is x(H 2 Oi ce ) ~ 4n gr <Tg r N s ^H 2 o/ n ~ 8 x 10~ 6 (cthV2 x 10~ 21 
cm 2 ). At this point, the water ice abundance increases rapidly with depth. Here, the 
photodesorption rate decreases with depth due to attenuation of the FUV flux by grains, 
but the rate of freezeout only decreases when the gas phase O abundance decreases. For the 
gas phase O abundance to decrease significantly requires many monolayers of ice. This is 
seen in Fig. 3 as the increase in x(H 2 Oj ce ) by two orders of magnitude at Ay ~ 3. Note the 
agreement of this freezeout depth with Ayj given in Eq. (16). As the gas phase abundance of 
O decreases, the growth of water ice slows and the ice abundance saturates (i.e., all O nuclei 
locked in water ice) near Ay ~ 3.5. Fed by photodesorption, the gas phase abundance of H 2 
maintains a broad plateau x(H 2 0)~ 10~ 7 from Ay ~ Ayf ~ 3 to Ay ~ 6, beyond which (i.e., 
deeper in the cloud) other gas phase H 2 destruction mechanisms besides photodissociation 
become dominant. The abundance of gas-phase H 2 near the peak is determined by the 
balance of photodesorption from the surface layer of ice on the grains with photodissociation 
of H 2 by FUV photons. Both of these processes are proportional to the local attenuated 
FUV field; this leads to a plateau in the H 2 abundance independent of Go (see §3.3). 

The abundances of other O-bearing molecules closely track the abundance of H 2 0. 
OH is mainly a product of photodissociation of gas-phase water, while it is destroyed by 
photodissociation and by reaction with S + , both of which depend on the attenuated FUV 
field. 2 is produced primarily through the reaction O + OH — > 2 + H and is destroyed 
by photodissociation. Thus the gas-phase H 2 0, OH and 2 abundances track one another 
from the beginning of freezeout through the H 2 plateau. We discuss this in more detail in 
§3.3 below. 

The Freezeout Layer. Beyond Ay > 6, photodesorption of water ice ceases to produce 
enough gas phase H 2 or OH to maintain the H 2 and OH plateaus (and therefore the 
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O2 plateau). Other destruction mechanisms than photodissociation of H2O take over, and 
gas phase abundances of all O-bearing species drop as H 2 ice incorporates essentially all 
available elemental 0. 

The destruction of CO by He + , followed by freezeout of the liberated O atom, as de- 
scribed in §2.3, is responsible for the rapid drop in the steady-state gas phase CO abundance 
at Ay ~ 6 (see Fig. 4) as all of the oxygen becomes locked in water ice. In our rather 
limited steady state carbon chemistry, C atoms end up incorporated into CH 4ice once O is 
liberated from CO. In our time-dependent runs, which have more extensive chemistry, the 
standard case at Ay = 10 has ~ 0.6 of the elemental carbon in CH 4 j ce and ~ 0.4 in COj ce 
after 2 x 10 7 yearsjf The ratio of the CH 4 j ce to the H20j ce abundance is ~ 0.15 — 0.25 at this 
time. Oberg et al (2008) discuss recent observational results on the ice content of grains in 
molecular clouds with emphasis on CH 4ice ; they conclude that CH 4ice does indeed form by 
hydrogenation on grain surfaces, as we have assumed in our models. However, they observe 
abundance ratios of CH 4 j ce to H20j ce that are typically 0.05, nearly 3-5 times smaller than in 
our standard model. This likely arises because our time dependent code starts with carbon 
in atomic or singly ionized form, so that there is more opportunity for C or C + to transform 
to CH 4 j ce than if the carbon started in the form of CO. Aikawa et al (2005) provide theo- 
retical models that result in the observed ratios of CH 4 j ce /H 2 0j ce by initiating carbon in the 
form of CO. The purpose of this paper is to follow H2O and 2 chemistry, and not carbon 
chemistry. We have therefore also neglected grain surface reactions that could lead to C02« C e 
or methanol ice@ Our C02i ce is formed by first forming gas phase CO2 by the reaction of 
gas phase CO with OH, and then subsequently freezing the CO2 onto the grain surfaces. We 
plan to study the possible fractionation signature ( 12 C02i ce / 13 C02ice) of our gas phase route 
to forming C02i ce in our subsequent paper which examines carbon chemistry in more detail 
(in preparation). In summary, our carbon chemistry here is simplified, and the abundances 
of carbon ices in the models are very crude. 



3.2. Dependence on Go, n, Yh 2 Oi °Hi PAHs, and b\ 

Dependence on Incident FUV Flux, Go- Fig. [5] shows the effect of changing the FUV 



5 The drop in the total CO gas + COj ce abundance by a factor ~ 2 in the opaque interior compared 
to the intermediate shielded regions would be hard to measure without a precise (better than factor of 2) 
independent measure of the column along various lines of sight through the cloud. 

6 We note that under optimum conditions C02i ce and methanol ice may carry nearly as much oxygen as 
water ice. However, generally they do not; therefore, their neglect should not appreciably affect our model 
of basic oxygen chemistry. 
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field strength, Go, on the depth of the gas-phase water peak and the abundance at the peak, 
while keeping the density fixed. Here Go is varied from 1 to 10 3 . The main effect of changing 
Go is to move the position of the water peak inwards with increasing Gq. Varying Go from 1 
to 10 3 changes the depth of freezeout from Ay ~ 1 to Ay ~ 7. Gas phase H 2 peaks once ice 
has formed at least a monolayer on the grain surfaces, so the location of the water peak also 
moves inwards (logarithmically) with increasing Go- The peak gas-phase water abundance 
(~ 1CT 7 ) is insensitive to G and insensitive to depth (i.e., A v ) in the water plateau. We 
also see that the depth where the gas phase water starts to decline from its peak (plateau) 
value is also weakly (logarithmically) dependent on Go, increasing with increasing Go- Thus, 
the column of gas phase water remains fairly constant and independent of G . 

However, some differences arise in the case with G = 10 3 . Here the assumption that 
every O sticking to a grain forms water ice breaks down. The high FUV field absorbed at 
the cloud surface leads to a high IR field that heats the grains to T gr ~ 25 K even at Ay ~ 8, 
so that a significant fraction of O atoms are thermally desorbed from grains faster than they 
can form water ice, resulting in a higher gas phase abundance of atomic O. At these large 
depths, the formation of gas phase water is not mainly by photodesorption of water ice, but 
by cosmic ray initiated ion-molecule reactions (see Fig. 1). The high gas phase elemental 
abundance of O leads to a high rate of formation of H 2 and 2 , whereas the high Ay leads 
to low photodestruction rates of these two molecules. In addition, although the FUV field 
is severely attenuated, the residual FUV is sufficient to prevent the complete freezeout of 
gas phase H 2 0. The combined effect is a deeper and higher peak in the gas phase H 2 and 
2 abundances. At very high Ay {Ay > 8) the freezeout of water ice finally drains the gas 
phase oxygen, and gas phase H 2 and 2 decrease with increasing Ay. Our time dependent 
results at t — 2 x 10 7 years find most of the oxygen in H 2 0-ice (~ 70%) and C0 2 -ice (~ 30%) 
at Ay > 8, and the carbon is mostly in CsH-ice (~ 60%) and C0 2 -ice (~ 40%) since the 
warm grain temperatures lead to thermal desorption of CO-ice and CH 4 -ice. 

Dependence on Gas Density, n. In Fig. 6 we vary the gas density from 10 3 to 10 5 cm~ 3 , 
keeping the rest of the parameters fixed at their standard values. Comparing this figure 
with Fig. 5 we see that for Go < 10 3 the plateau water abundance is independent of 
gas density and G , depending only on the photoelectric yield Y H , 2Q and the grain cross 
sectional area per H nucleus a H as shown below. The peak abundance is independent of 
n because both the formation rate per unit volume (photodesorption) and the destruction 
rate (photodissociation) are proportional to n. For a fixed value of Go, lowering the gas 
density moves the water peak further in; the deposition of water ice depends on the product 
n gr no oc n 2 , while the photodesorption rate goes only as the grain density n gr oc n times the 
local FUV field. Therefore, the water peak moves into the cloud with decreasing density. 
However, the column or Ay where the water starts to decline from its peak (plateau) value 
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is insensitive to n. Thus, the plateau starts to become narrower and narrower (in AAy) as 
n decreases. As a result, The column of gas phase water increases with n somewhat. We 
explain all these trends in more detail below in §3.3, where we present a simple analytic 
analysis. 

Dependence of Gas Temperature on Go and n. Figure 7 presents the gas temperature 
as a function of Ay for a range of our Go, n parameter space. This figure reveals the 
temperature at which most of the H 2 and 2 is radiating, as we have marked the peak in 
the gas phase water abundance (which coincides with the peak in the O2 abundance) as solid 
portions of these curves. We see that over the parameter space Go — 1 — 10 5 the emission 
from H2O and O2 arises from regions of T ~ 7 — 50 K, with T generally rising gradually 
with G . For G < 500, there is a true "water plateau" and the solid line corresponds to 
the plateau region. T ~ 7 — 30 K in the plateau, and the gas temperature is often higher 
than the dust temperature because of grain photoelectric heating of the gas. For higher 
Go there is no true plateau but the gas phase H2O and O2 peak at higher Ay, where the 
gas and dust temperatures are nearly the same. We will show in §3.4 that at high Ay the 
dust temperature scales roughly as G[}' 24 . Here we have made solid lines where the H2O 
abundance is more than 0.33 of the peak abundance. 

Dependence on Photodesorption Yield, Yh 2 o- Fig. 8 presents the profiles of the gas phase 
H 2 and O2 abundances, and the H 2 ice abundance as the photodesorption yield Yh 2 o of 
water ice is varied, maintaining Y h, w = 2Yff 2 o, but keeping all other parameters standard. 
In contrast to the lack of dependence of the plateau value of the gas phase water abundance 
x p i(H 2 0) on n and Go, x p i(H 2 0) scales roughly linearly with Yh 2 o for low Yh 2 o < 10~ 3 . 
This linear dependence can be easily understood: the formation rate (photodesorption) of 
gas phase H 2 scales linearly with Yp 2 o whereas the destruction rate (photodissociation) at 
a given Ay is independent of Yff 2 ol3 On the other hand, at low values of Yh 2 o < 10~ 3 , the 
abundance of 2 rises with Y£ Q , or, in other words, with [x(H 2 0)} 2 . 

Dependence on Effective Grain Area, an- In Figure 9 we show the effect of varying 
on with all other parameters at their standard values. The abundance of gas phase H 2 
in the plateau x p i(H 2 0) increases nearly linearly with as- The depth Ay ~ Ayf where 
the H 2 abundance first reaches its plateau value is insensitive to cr#, as predicted by Eq. 
(17). Similarly, the depth where H2O finally begins to drop from its plateau value also is 



7 Where thermal desorption of atomic O can compete with photodesorption, i.e. high grain temperatures 
from high FUV fields, this is not the case. In addition, at high values of Yh 2 o > 10~ 3 , the surface can 
become depleted in water ice compared with, for example, CO ice, which then lowers the desorption rate 
of H2O into the gas. This leads to a smaller increase in gas phase H2O with increasing Yh 2 o than would 
otherwise be the case. 
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insensitive to er#. The abundance of O2, like the abundance of H 2 0, scales with <r#. 

Dependence on PAHs. We are most interested in the chemistry at Ay ~ 3 — 6 where 
the H 2 and 2 abundances peak. Because it is not clear if PAHs exist at such depths 
(Boulanger et al 1990, Bernard et al 1993, Ristorcelli et al 2003), we have not included 
PAHs in our standard case but have assumed a grain surface area per H given by on which 
is appropriate to an MRN distribution extending down to 20 A grains when computing the 
neutralization of the gas by collisions with small dust particles. If PAHs are present, the 
effect of neutralization is enhanced, which changes the gas phase chemistry somewhat. We 
have found that the inclusion of PAHs into our standard case does not affect the results 
appreciably. There is a somewhat greater effect at lower density such as n ~ 100 cm -3 , 
where the inclusion of PAHs does increase the total water column to high Ay by a factor of 
about 10, both by broadening the peak plateau and by increasing the peak abundance (for 
the case Go = !)• 

Dependence on b\. We have run our standard case but have assumed that the small 
grains responsible for FUV extinction are depleted, for example, by coagulation. In our 
test run we have decreased all our b x by a factor of 2, reduced cr# by 2, and reduced grain 
photoelectric heating accordingly. The main effect is to increase the characteristic Ay where 
the H2O and O2 plateau by a factor of 2. However, the abundances of H2O and O2 in the 
plateau decrease because of the decrease in an- Therefore, the columns of gas phase H 2 and 
2 only increase by a factor of about 1.5 from the standard case. The gas temperature does 
not appreciably change. In summary, the change in the intensities and average abundances 
of H 2 and O2 is not large, and is expected to be less than a factor of 2. 

Dependence of H2O and O2 Columns on Gq andn. In Fig. [10], we show the total columns 
of H 2 and O2 for a cloud with high Ay > 10 as a function of the incident FUV flux Go and 
the cloud density n. For low G < 500, the H 2 and 2 columns are roughly independent of 
Go and n, and are of order 10 15 cm -2 for our standard values of photodesorption yield. These 
columns arise in the H2O and O2 plateaus, at intermediate cloud columns. In the steady state 
models, there is so little gas phase water in the interior of the cloud that the interior does 
not contribute significantly to the total column. Therefore, once the total column through 
the cloud iV is large enough to include the plateau (N H > 10 22 cm -2 ), the average abundance 
will scale as N# _1 . Thus, part of the reason for the low average abundances observed by 
SWAS and Odin is the dilution caused by the lack of H 2 and O2 at the surface and deep 
in the interior of very opaque, high column, clouds. At higher G > 500 we see the effect of 
the thermal evaporation of the O atoms from the grain surfaces. The effect is most dramatic 
for the gas phase O2 column, which rises to values N(C>2) ~2 x 10 16 cm -2 . This is close to 
what is needed for a detection by Herschel and therefore it is in this type of environment 
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that Herschel may potentially detect O2. We note, however, that this result depends on the 
binding energy of O atoms to water ice surfaces. The value we have adopted (see Table 1) 
is standard in the literature but it refers to van der Waal binding to a chemically saturated 
surface. It is quite possible that the O binding energy is larger, and this would then increase 
the temperature (and Go) where this effect would be initiated. 



3.3. Simple Analytic Analysis of the Results 

The results in §3.1 and §3.2 can be understood by a simple analytic chemical model 
that incorporates the main physics. Such a model, though approximate, has the advantage 
of allowing one to determine and understand the sensitivity to various model parameters, 
and serves to validate the numerical model. In this subsection the chemical rates are all 
taken from the UMIST database, except for the ion dipole rates and the photodesorption 
rates discussed in §2.7. In the following, the reader is referred to Table 3 in the Appendix 
for quick reference to symbols. 

The abundance x(H20) of gas phase H2O can be roughly traced and understood by 
simple analytic formulae. Photodissociation dominates the destruction of gas phase H 2 
from the surface to the depth A Vd where gas phase destruction with B.^ generally takes over. 
From < A v < A Vd the simple model assumes that H 2 is formed by photodesorption of 
H 2 ice. Therefore, by equating formation to destruction 

GvF Q e~ l * Av Y H2 ofs,H 2 on 9r (J g r = G Q R H20 e~ l - 1Av n(H 2 0), (19) 

or 

x(H 2 0) = FoY^oUHzOVHRHloe' - 1 ^, (20) 

where Rh 2 o — 5.9 x 10~ 10 s -1 is the unshielded photodissociation rate of H 2 in an FUV field 
of G = 1. The differences in the b\ factors of 1.7 and 1.8 is because the two processes are 
dominated by somewhat different FUV wavelengths. The fractional coverage of the surface 
by water ice f s ,H 2 o is given by equating the sticking of O atoms to grain surfaces to the 
photodesorption rate of H 2 (including both desorption of H 2 and of OH + H). Therefore, 
recalling that Y = Y H2 o + Y OH , w , 

n(0)v n gr a gr = G F Q e~ 18Av Y f s ,H 2 on gr (T g r (21) 
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or 

_ n(Q)v e l - SA v 
Uh2 ° ~ G F Y (22) 
Once a monolayer forms, f s ,H 2 o ~ 1- We have found in our numerical runs that sometimes 
fs,H 2 o saturates at ~ 0.5 rather than 1.0, because of the presence of other ices such as CO 
ice or CH 4 ice. However, assuming a saturation value of unity, the critical depth A V f at 
which a monolayer forms is then given (this is equivalent to Eq. 17 in §2.5): 

A Vf = 0.56\n(^f-), (23) 



n(0)v 

or a critical H nucleus column (Spitzer 1978) from the surface of 

Nf = 2 x 10 21 A Vf cm- 2 . (24) 
We therefore find, using equations (20 and 22), 

x(H 2 0) = — — for A v < A vf . (25) 

I Lr ItH 2 O 



Xpl {H.O) = toiJ^l . 2x10- ^ ., , 'I ., io,A Vf < A v < A Vd . 



F Y H 2O a H e-^ ^ (Y H2 p\ ( a_n 

R H2 o ~ V10- 3 ; V 2 x 10- 21 cm 2 

(26) 

The abundance of gas phase H2O rises exponentially with depth (as e L7j4v ) until a monolayer 
of water ice has formed on the grains. It then forms a plateau with value x p i(H 2 0) that is 
independent of Go and n and weakly dependent on Ay. The plateau abundance goes linearly 
with Y and an- The lack of dependence on Go can be understood since both the formation 
and destruction of gas phase H 2 depend linearly on Go- Higher fluxes do decrease the 
surface abundance of gas phase water, but the freezeout depth is deeper, so there is more 
dust attenuation, and the gas phase water abundance rises and peaks at the same constant 
value. More precisely, regardless of the incident FUV flux, the local attenuated FUV flux at 
the freezeout point is always the same, because it must be the flux which photodesorbs an ice 
covered grain at the rate that oxygen from the gas resupplies the surface with ice. Therefore, 
the gas phase water abundance is the same as well, regardless of Go, only depending on 
the the assumed values of the grain parameters and the photodesorption yield. The lack 
of dependence on n can be understood similarly since both the formation rate per unit 
volume and the destruction rate per unit volume depend linearly on n. On the other hand, 
x(H 2 0) oc Ycth, since the formation rate by photodesorption goes as this product whereas 
the destruction rate is independent of them. 
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In order to determine Aya and to provide an approximation for x(H 2 0) that applies 
for Ay > A Vd , we need to include other destruction routes for gas phase H 2 that begin to 
dominate at high Ay. We find that generally the route which takes over from photodissoci- 
ation is reaction with Hjj" to form H 3 + . H 3 + recombines one fourth of the time to reform 
water (not a net destruction) and three fourths of the time to form OH (Jensen et al 2000, 
Neufeld et al 2002). The net rate coefficient for the destruction of gas phase H 2 by H3 is 
then 



1h+ = 0.75(4.5 x 10 _9 T. 



-1^ 



300 



cm 3 s 
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where T300 = T/300 K. Equating the formation of gas phase H2O by photodesorption to the 
destruction of H 2 by both photodissociation and Hjj", we obtain a more general solution for 
x(H 2 0) for Ay > Ayf that applies beyond A Vd 



x(H 2 0) 



x P i(H 2 0) 



1 + 



3 

GoRh 9 o 



for Ay > Ay f . 



(28) 



Note that the second term in the denominator, which represents destruction by hLt , begins 
to dominate when its value reaches and exceeds unity. At this point, x(H 2 0) begins to drop 
from its plateau value. 

Hjj" is formed by cosmic rays ionizing H2 followed by rapid reaction of with H2. 
When gas phase CO is depleted, it is destroyed by dissociative recombination with electrons. 
The electrons are provided by Si + and S + at moderate Ay and by H3 , He + , and H + at high 
A v , so the electron density is complicated. At high Ay the code results suggest that the 
electron density is roughly n(e) ~ 2n(if 3 l ~). We then obtain 



n(H+) highAv ~ lO-^V/ 2 cm" 3 , (29) 

with H nucleus density n in cm -3 . In Eq. (29) we have assumed a total (including secondary 
ionizations) cosmic ray ionization rate of 5xl0~ 17 s -1 for H 2 (Dalgarno 2006). At moderate 
A v Si + usually provides most of the electrons. Gas phase Si + is produced by photoioniza- 
tion of gas phase Si and destroyed by recombination with electrons. The Si + density, and 
therefore the electrons provided by Si + , exponentially drops with Ay as the FUV photons 
that photoionize Si drop with increasing dust attenuation. As long as most of the gas phase 
silicon is Si + , the electron density is high and this suppresses the abundance of hht , prolong- 
ing the water plateau to higher Ay (see Eq. 28). However, we find that once the gas phase 
Si + density drops to 2n(H^) high A v , so that the electrons provided by Si + are comparable to 
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those provided by H3 , He + , and H + , n(H^) has risen to ~ n(H^) high A v and the H3 term in 
the denominator of Eq. (28) dominates. Here the gas phase water abundance starts to fall. 
This occurs when 



A Vd = 0.59 In 



4.5 x 10 5 G n- 1 / 2 T 30 ° - 2 (- 

V 



n Si T ' 



(30) 



UH 2 0-gr , 

where nsir = 1-7 x 10~ 6 n is the density of silicon both in the gas phase and as ice mantles 
on grains, and n^o-gr ~ 3 x 10~ 4 n is the density of water molecules incorporated in ice 
mantles. 

We have checked our analytic expressions for A V f (Eq. 23), A Vd (Eq. 30), and x pl (H 2 0) 
(Eq. 26) and find good agreement (Ay/ and Ay d better than ±25% and x p i(H 2 0) good to 
a factor of 2 over our (Go,n) parameter space for Go < 500). For Go > 500 the grains are so 
warm that O atoms adsorbed to grains thermally evaporate before reacting with H atoms to 
form OH on the grain surface. As discussed in §3.2, this results in a change in the behaviour 
of gas phase H2O. 

The abundance of O, OH, and O2 in the H2O plateau region Ayf < Ay < Ay d can also 
be analytically estimated from the above expression for x p i(H 2 0). We first find x p i(OH) by 
equating formation by photodissociation of H 2 and photodesorption of OH from water ice 
with destruction by photodissociation of OH and reaction with S + and He + . The photodes- 
orption rate of OH is 2 times the rate of H 2 photodesorption, and therefore 2 times the 
rate of H 2 photodissociation. On the other hand, the destruction of OH by S + and He+ is 
roughly 2 times that OH photodissociation. Therefore, 



x pl {OH) ~ [~^J x pi {H 2 0) ~ 1.7x pl (H 2 0), A Vf < A v < A Vd (31) 

where Ron = 3.5 x 10 -10 s _1 is the unshielded photodissociation rate of OH when Go = 1. 
Thus, OH tracks H 2 with somewhat greater abundance in the plateau, as observed in 
Figure 3. 

The gas phase O abundance is determined by balancing the formation of O by photodis- 
sociation of OH with the removal of O by adsorption onto grains: 



x p i{0) 



GoRoh& 



-1.7/H 



<j H nv 



x pl (OH) ~4 x 10" 



G /n 

10- 2 



,-1.7A v 
IO- 3 



x P i{H 2 oy 



10- 



(32) 



Therefore, x p i(0) exponentially declines with e 1 - 7Av due to the declining photodissociation 
of OH, but the constant (with A v ) collision rate of O with grains, as seen in Figure 3. It 
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also increases with Go/n for fixed Ay. This formula breaks down when the RHS exceeds 
~ 10~ 4 , where x p i(0) saturates since all elemental gas phase O is now in atomic form. 

The O2 abundance follows from the abundances of O and OH since O2 is formed by 
the neutral reaction of OH with O (rate coefficient 7o 2 = 4.3 x Kr n (300 K/T) 1 / 2 e" 30 K/T 
cm 3 s _1 ). O2, like H2O and OH, is destroyed by photodissociation in the plateau region. 
Therefore, 



x p i{0 2 



l-7 2 7o 2 
ohVq 



V K o 2 J 



e°- 1A -[x pl (H 2 0)] 2 (33) 



2 x 10 21 cm 2 \ /30 x r ' " ^- l2 
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where Ro 2 — 6.9 x 10 -10 s _1 is the unshielded rate of photodissociation of O2 in a Go = 1 
field. The [x p i(H 2 0)} 2 dependence, which translates to a dependence in the plateau for 
low Yh 2 o < 10~ 3 , is clearly seen in Figure 8. Note also that Eq. (33) explains the lack of 
dependence of x p i{0 2 ) on n and G , and also predicts that x p i{0 2 ) oc a H as observed in the 
numerical models. 

In summary of the above, the approximate analytic equations can be used in the cases 
with Go ^,500 to understand and quantitatively estimate: (i) x p i(H 2 0,OH,0 and O2) and 
their dependence (or independence) on n, Go, Y, and an] (ii) the depth Ayj of the onset of 
the plateau, the depth A Vd of the termination of the plateau, the width AA V = A Vd — A V f of 
the plateau, and their dependence on the above parameters. At high Y H2Q > 10~ 3 the abun- 
dances of O, OH, O2, and H2O are overestimated by factors of order 2 in the above formulae 
because we have ignored the presence of other ices, which reduce the fractional surface cov- 
erage of water ice, and therefore the production of gas phase H 2 and its photodissociation 
products. 

Finally, the above analytic expressions can be used to predict the optimum conditions 
for observing H2O and O2. In general, for our model runs the intensity of the H2O ground 
state ortho line at 557 GHz and para line at 1113 GHz can be estimated by assuming that the 
emission is "effectively thin" Jf| In other words, although the optical depth in the ground state 
may be larger than unity, the gas density is so far below the critical density (~ 10 8 cm -3 ) 
that self absorbed photons are re-emitted and eventually escape the cloud (see discussion in 
§5). In this limit every collisional excitation results in an escaping photon, and 



8 We include the 1113 GHz line because it is accessible with the upcoming Herschel Observatory. 
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hv[x nj eXi0 + x p nj ex , p ]N(o - H 2 0) 
^ ^- , (34) 

where 7 eXj0 ~ 2.5 x \Q-^ e ~ 2& - 7 / T and 7 ex>p ~ 3.5 x \Q-n e -2&.7/T cm 3 g -i ( valid for t<30K) 
are the rate coefficients for the collisional excitation of the 557 GHz transition by impact with 
ortho and para H 2 respectively (Faure et al 2007; M-L Dubernet, private communication). 
Here, x (x p ) is the abundance of ortho(para)-H 2 with respect to H nuclei, hu is the energy 
of the 557 GHz photon, and N(o — H 2 0) is the column of ortho water. Similarly, the para 
H 2 1113 GHz line intensity is 

hv[x n~f eX)0 + x p nj ex>p )N(p - H 2 0) 

J 1113 — J- , [oO) 

where where ^ exfi ~ 3 X ]_0-lOg--53.4/T ^ g q x ^g-ll e -53.4/T cm 3 s -l (ygJid f or 

T< 30 K) are the rate coefficients for the collisional excitation of the 1113 GHz transition 
by impact with ortho and para H 2 . From the above equations, assuming most of the H 2 
column arises in the plateau between Ayf and Ayd, the ortho(para) H 2 column is given 

N(o(p) - H 2 0) ~ 2 x ltf\A Vd - A Vf )x pl (H 2 0)f (f p ) cm" 2 , (36) 

where f (f p ) is the fraction of H 2 in the ortho(para) state. Substituting our analytic 
equations for these parameters, we obtain 

^ »*Hr T /.Mo.n*] (w^) {wS) ^-°-""e— /. erg cm - 2 s -. sr -> 

(37) 

and 

- 2A (!) - 26,/Tfe " < 38 > 

where Av4y = Ayd — A v f. Note that the column of H 2 is relatively independent of Go and n. 
However, the intensity I 557 or I lll3 is proportional to the H 2 column times ^ e x n oc ne~ 26 ' 7 ^ T 
or ne _53 ' 4//T , and therefore is somewhat sensitive to Go (which determines T) but is more 
sensitive to n. The 557 GHz line depends on f Q which is a sensitive function of T gr for low 
Go (i.e., low T gr ), as will be described below. Therefore, the 557 or 1113 GHz H 2 line will 
be strongest in dense regions illuminated by strong FUV fluxes (assuming sufficient column 
or Ay to encompass the plateau, and assuming the regions fill the beam). We note that the 
gas and grain temperatures vary somewhat through the water plateau region. We take the 
gas T and the grain T gr to be the average of the values at the freezeout point Ayf and at 



Note that the maximum values of x Q and x v are 0.5. 
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the Ay where the gas phase water abundance peaks. Figure 11 shows a contour plot of this 
average gas temperature T ave as a function of Go and n. Over nearly the entire range of Go 
and n, T ave only varies from ~ 15 — 40 K. Recall that Eqs (37) and (38) break down for 
Go ^ 500, but note that from our numerical results -Z557 and I\uz will increase further, and 
the analytic equations will underestimate the intensities with increasing Go above Go ~ 500, 
since the water column rises here. 

Similarly, we can estimate the intensity of the O2 transition. Here, the density is above 
the critical density (~ 10 3 cm~ 3 ) for the 487 GHz O2 transition so the levels are in LTE. In 
addition, because the transition has a low Einstein A value, large columns are needed to make 
the transition optically thick. Thus, again assuming optically thin emission, the intensity 
J 487 oc N (0 2 )e~ 25K / T / Z (T) , where Z(T) is the partition function. An approximation to 
Z(T) for K< T < 25 K is given as Z(T) ~ 3 + 2.1T + 0.1T 2 . Using the equations from 
this subsection, we then estimate 



J 487 ~ 6xHT 8 ( ^] ( (^4) ( %-^] e-^ T e-°-°^f +A ^ erg cm" 2 s" 1 sr" 

87 \Z{T)) V T 7 V !0 3 / V2xl0- 21 cm 2 y & 

(39) 

Note that again, this equation breaks down for G ^ 500, and the equation will underestimate 
the intensity for higher G since the 2 column rises. In general, however, the equation 
predicts little n dependence (there is a small n dependence in the Ay terms). At low 
Go ^ 500, the equation predicts only moderate Go dependence. Although various terms 
depend on T, the temperature of the (3 2 plateau region is not sensitive to Go since the plateau 
retreats further into the shielded region for higher Go- In addition, the T dependence in the 
exp(-55/T) term is partially cancelled by the TZ(T) term in the denominator. 



3.4. The H 2 557 and 1113 GHz and 2 487 GHz Line Intensities 

Using the numerical models, Figures 12 and 13 present H 2 -Z557 and /1113 and 2 1437 
contours as functions of n and Go for fixed standard values of the other parameters. We 
have assumed that the H 2 ortho to para ratios are in LTE with the grain temperature 
(which in the plateau is nearly the same as the gas temperature except for very low ratios 
of Go/n, which bring the plateau very close to the surface of the cloud where T > T gr ). The 
grain temperature is appropriate since a given H 2 molecule forms on a grain and spends 
too little time in the gas to equilibrate to the gas temperature before it is photodissociated. 
The grain temperature in the plateau is approximately (Hollenbach et al 1991) 
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~ (3 x 10 4 + 2 x 10 3 G 4 ' 2 ) 



(40) 



A fit to the LTE ortho to para ratio of H 2 valid to 10% in the region 5 K < T gr < 20 K 
and good to 20% for T gr > 20 K (e.g. Mumma et al 1987) is given: 



The H2 ortho to para ratio is calculated in steady state using the rates and processes 
described by Burton et al (1992). They are neither 3 nor LTE. They are insensitive to n in 
the plateau, but sensitive to Go, which affects T. The H 2 ortho to para ratio x a /x p ~ 0.5 
at G = 3 x 10 4 , 0.3 at G = 1.5 x 10 4 , and 0.1 for G = 4x 10 3 . For lower values of G Q , 
x jx p < 0.1 so that para H2 collisions dominate the excitation of water. 

The 557 GHz and 1113 GHz intensities behave as described by the analytic formulae 
(Eq. 37 and 38) presented in §3.3 for Go ^ 500. The H 2 intensities rise roughly linearly 
with n. For fixed density ^557 and /1113 tend to rise with Go (which raises T gr as shown in 
Eq.(40) and which tends to raise T as can be seen in Fig. 11) partly due the excitation 
factors e~ 26 7//T and e _53 4//T , but also because increasing T increases the ortho/para ratio of 
H2, and the ortho excitation rates in collisions with H2O are much higher. /1113 does not rise 
as much with increasing Go as one might naively expect given its excitation factor e _53 - 4 / T , 
because the increasing T reduces f p from a value of ~ 1 at low G to 0.25 at high G o 

The following summarizes the accuracy of the analytic approximations for the H 2 and 
2 intensities for Y H , z0 ~ 10~ 3 . For G < 500 the analytic expression for I 557 agrees with 
the numerical results to within a factor of 3 for n = 10 3 — 10 5 cm -3 . The agreement for 
/1113 is within a factor of 3 for n = 10 4 — 10 5 cm -3 , but only to a factor of 5 (with the 
analytic underpredicting the intensity) for n = 10 3 cm -3 . However, for this low value of n, 
the intensities are weak and not detectable. For higher values of Go the analytic formula 
underpredicts the intensity because of the larger columns of H 2 and 2 (see Figure 10 
and discussion). For higher values of n the water plateau is so close to the surface that the 
gas temperature drops rapidly and significantly from Ayj to Ayd] here the analytic formula 
is not accurate since it assumes constant average temperatures for both grains and gas. 



10 Note also that because increasing Go pushes the water and O2 plateaus to higher Ay, the gas temperature 
T ave only rises weakly with Go However, J557 increases with increasing / // p , which is very sensitive to 
increasing T gr (or Go) for low T gr (or Go). O2 behaves somewhat differently, being less sensitive to n, as 
predicted by the analytic models. 
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Nevertheless, the analytic formulae are useful in making estimates of line intensities, and 
also illustrate the dependences of the intensities on the physical parameters. 

3.5. Model of Chemical Desorption 

Following O'Neill, Viti, & Williams (2002), as discussed in §2.4, we have run our stan- 
dard case with the modification that all molecules formed on grain surfaces instantly desorb. 
In our standard case only H 2 instantly desorbs. The basic differences are that the freezeout 
depth Ayf moves inward (increases) by about 2 magnitudes of visual extinction due to the 
chemical release of adsorbed OH and H 2 0, the peak gas phase abundance of H 2 increases 
by about a factor of 10, and that of 2 by a factor of roughly 100. We do not pursue this type 
of model further here, but note that it could warrant further study, although it appears to 
predict excessive H 2 and 2 intensities. In the original O'Neill et al work, photodesorption 
and cosmic ray desorption were not included. 

3.6. Model with No Grain Chemistry 

In order to test the importance of the formation of OH, H 2 and carbon molecules on 
grain surfaces in our model, we have run the standard case (Go — 100, n = 10 4 ) and two 
other cases (G = 100, n = 10 3 cm" 3 and G = 10, n = 10 4 cm" 3 ) with only H 2 formed 
on grain surfaces and no other grain chemistry. Like the model described above in §3.5, 
one main effect is to push the peak in the H 2 and 2 abundances to higher Ay, since 
grain formation of H 2 is suppressed. However, in this case there is no enhancement of the 
peak abundances, because the formation rates of OH and H 2 are suppressed even more 
than in the case studied in §3.5. In addition, there can be significant depletion of gas phase 
elemental O at low Go by the formation of atomic O ice, which does not occur in §3.5 or in 
our models with grain chemistry in §3.1 and §3.2. In the (Go = 100, n = 10 4 cm" 3 ) case 
and the (Go — 10, n — 10 4 cm" 3 ) case, there is a reduction in both the columns of H 2 and 
2 and in the gas and grain temperatures of the emitting region compared to the cases with 
grain chemistry. Thus, the predicted lines are significantly weaker. Because of the good 
quality of our grain chemistry model fits to the data described in §4, and because of the 
laboratory and observational evidence for grain chemistry described in §2.3 and §3.1, we do 
not consider this case further. 
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3.7. Time Dependence 

The numerical results and analytic solutions presented so far are for steady state solu- 
tions to the chemistry. However, as noted in §2.3 and §2.4, in the opaque interior of molecular 
clouds at A v > A V d, the chemical timescales may exceed or be comparable to the cloud age, 
so that time dependent solutions to the chemistry are needed. The long chemical timescales 
include the freezeout of gas species onto grains for n<10 3 cm" 3 (see Eq. 1), the cosmic 
ray desorption of CO ice (see §2.4), and the conversion of gas phase CO to water ice which 
is initiated by He + (see §2.3). In the opaque "freezeout" regions the attenuated FUV flux 
plays little role in the chemistry so that the time dependent abundances do not vary with 
Ay. Therefore, we present results at a representative Ay = 10 for the abundances of gas 
and ice species as a function of the age of the cloud, the density of the cloud, and the cosmic 
ray desorption rate of CO. 

The time dependent code is that of Bergin et al (2000), updated to have the same grain 
chemistry and desorption processes as described above for the steady state code. It includes 
more species than the PDR surface code, and therefore tracks the carbon chemistry more 
accurately in the inner regions. It takes as initial conditions the chemistry of translucent 
clouds: all the hydrogen is initially molecular H 2 with the other species atomic (for ionization 
potentials > 13.6 eV) or singly ionized otherwise. 

The results are shown in Figures 14 and 15, where we present the abundances of O2, 
H 2 and CO plotted as a function of the cosmic ray desorption rate of CO ice, for two times 
representative of cloud ages (2 x 10 6 years in Fig. 14 and 2 x 10 7 years in Fig. 15), for 
the standard G = 100, and for 3 densities (10 3 , 10 4 , and 10 5 cm -3 ). Of particular note is 
Figure 15, where one sees the gas phase molecules increase in abundance as the cosmic ray 
desorption rate increases, reach a peak when the timescale for cosmic ray desorption of CO 
is about 0.1 - 1.0 times the cloud evolution time of 2xl0 7 years, and then decreases as the 
cosmic ray rate increases further and the desorption time becomes less than ~ 2x 10 6 years. 
If the rate is too low, then the CO ice remains on the grain and the gas phases abundances 
are low. If the rate is too high, then the CO ice rapidly desorbs, but then the CO in the gas 
breaks down by reaction with Hc + and the gas phase elemental oxygen freezes out as water 
ice, thereby lowering the gas phase abundances of CO, H 2 0, and 2 . The main point is that 
in all cases for t = 2 x 10 7 years and for most cases with t = 2 x 10 6 years, the abundances of 
H 2 and 2 at high A v , even when time dependent chemistry is taken into account (steady 
state chemistry is more extreme in this respect as seen in Figure 3), are substantially lower 
than the peak abundances of these molecules at intermediate Ay. Thus, most of the column 
of these molecules in molecular clouds are not provided by the interior regions, where they 
are frozen out, nor by the surface, where they are photodissociated. Rather, most of the 
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column arises from the intermediate zone. 

3.8. Chemistry: the Gas Phase C/O Ratio 

Figures 16 and 17 show the gas phase elemental abundance of C and the gas phase 
abundance ratio of elemental C / 'O , as a function of Ay, for two different cloud ages, t = 2 x 10 6 
years and t — 2 x 10 7 years, respectively. In each figure, we plot the standard case, the 
standard case but with the density changed to n — 10 3 and 10 5 cm -3 , and the standard 
case but with the incident FUV field changed to Go = 1 arid 10 3 . For these results, we have 
used our time dependent models, since the interesting effect of C/O exceeding unity occurs 
at high Ay, where the time dependent models are more relevant and are more accurate in 
treating the carbon chemistry. At early times, t ~ 2 x 10 6 years, the C/O ratio generally 
stays below unity at high Ay. The exception is the high density (n = 10 5 cm -3 ) case, 
but here the gas phase elemental abundance of carbon is extremely low, which implies low 
abundances of any organic molecules produced in this carbon-rich environment. At late 
times, t ~ 2 x 10 7 years, we see that all cases except the high Go = 10 3 case result in 
gas phase C/0> 1 for Ay > 8, due to the preferential freezeout of water ice which robs the 
gas of elemental O. For those cases with C/0> 1, the gas phase elemental C abundance 
increases with increasing Go and decreasing n, due to increased thermal desorption rates of 
carbon molecular ices and decreased freezeout rates, respectively!"] Figure 17 shows that for 
densities n ~ 10 3 — 10 4 cm -3 and Go = 1 — 100, strong chemical effects might be observed 
deep in clouds, with high abundances of gas phase hydrocarbons and other carbon molecules 
rather than CO. Unfortunately, these results can only be viewed as suggestive, because the 
gas phase C depends on the desorption rate of CO, CH 4 and other carbon species frozen to 
grain surfaces, and these are uncertain. Nevertheless, it appears that regions with C/O > 1 
and high gas phase elemental C abundance can appear deep in the cloud. This is consistent 
with the chemistry of GMC cores such as Orion A (e.g., Bergin et al 1997) and regions 
in the Taurus Molecular Cloud (e.g., Pratap et al 1997). In addition, in these deep layers 
significant amounts of CH 4 ice may form, and its release to the gas stimulates the carbon 
chemistry. In summary, the freezeout of the oxygen into water ice can have profound effects 
on the carbon chemistry deep in the cloud and the ability of the cloud to produce carbon 
chain molecules or hydrocarbons. 



11 At t = 2 x 10 7 years the time-dependent code shows that CO and CH4 ice dominates the carbon ices 
for Go ^500, and CO2 and C3H ice dominates at higher Go, due to thermal desorption of the more weakly 
bound CO and CH 4 . 
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3.9. Limb Brightening and Surface to Volume Effects 

Because our model suggests a peak in the H 2 and O2 abundance at intermediate 
depths into a molecular cloud, one might expect clouds to be limb-brightened if viewed 
with beams smaller than the cloud. We make here two caveats to this prediction. First, 
GMCs are clumpy, and the dumpiness means that Ayj will occur at different physical 
depths (i.e., radii) in the clouds, depending on the particular radial ray through the cloud 
and where it intersects clumps. This will have the effect of smearing the limb-brightening. 
In addition, although the 557 and 1113 GHz H2O transitions are "effectively thin", they 
do suffer numerous scattering events since the optical depth in the lines are greater than 
unity. This scattering will reduce the limb-brightening effect for the H 2 transitions. The 
O2 does not suffer scattering, so it may be more promising. However, the sensitivity of 
Herschel may require large O2 columns > 10 16 cm -2 , and these may only be achieved for 
high Go. In such cases, the peak O2 abundance occurs deep in the cloud, Ay ~ 8, and 
makes limb brightening only possible to detect in very high column but extended nearby 
clouds illuminated by high fields. Nevertheless, Herschel observations of O2 in clouds may 
find regions of limb-brightening. 

One still would like to map the clouds, however, even if limb brightening is hard to 
detect. The regions of high Go may be fairly limited in spatial extent, and our models 
predict higher H 2 and especially O2 columns in these regions which Herschel may spatially 
resolve. In addition, the H2O and O2 plateaus extend for AAy ~ 3. For gas densities of 
n ~ 10 4 cm -3 this corresponds to ~ 0.2 pc, which at 500 pc corresponds to an angular size 
~ 1.5'. Therefore, for a PDR slab viewed edge-on, Herschel may resolve the plateau region. 

Even for clouds too small in angular extent to map, our model can be tested by cor- 
relation studies of the H 2 and O2 intensities with the intensities of rotational transitions 
of other molecules known to trace either the surface or the volume of the cloud. For exam- 
ple, 13 CO and C 18 are thought to generally^ trace the volume of the cloud, whereas PDR 
species such as CN trace the "surface" of the cloud (Fuente et al 1995, Rodriguez- Franco et al 
1998). We are investigating such observational correlations in Melnick et al (in preparation). 



12 In some cases CO is frozen on grains in very high Ay regions at cloud centers; see, for example Bergin 
& Tafalla (2007) for a review of the observations. 
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4. APPLICATIONS TO SPECIFIC SOURCES 

4.1. Diffuse Clouds 

The diffuse clouds that we model, based on those observed by SWAS towards W51 and 
W49 where H 2 557 GHz is seen in absorption, are characterized by G ~ 1 — 10, n ~ 100 



cm , and Ay ~ 1 — 5 [these are sometimes called "translucent clouds" (ISnow &: McCall 



20061 )]. Diffuse clouds are illuminated on all sides by the diffuse interstellar FUV field, 
whereas our models are those of ID slabs illuminated from one side. To roughly construct 
a model of a diffuse cloud of total thickness A vt , we take the results from A v = to 0.5 Ayt 
and assume a mirror image for the rest of the cloud from 0.5Ay t to Ay t . 

We compare a single model fit of a diffuse cloud with the 54 and 68 km s _1 yelocity com- 
pone nts observed by SWAS towards the background submillimeter source W49 (IPlume et al. 



200 



3) . For 



20041 ) as well as in the single absorption component towards W51 (INeufeld et al. 
all three of these components, both OH and H2O column densities have been measured, 
with N(H 2 0)~ 1 - 3 x 10 13 cm~ 2 and N(OH)/N(H 2 0) ~ 3-4. Our diffuse cloud model with 
n = 10 2 cm" 3 and G = 1 gives N(H 2 0)=3 x 10 13 cm" 2 and N(OH)/N(H 2 0) =3.5 for a total 
cloud extinction of Ay t = 4.4 magnitudes. We find that assuming a photodesorption yield 
Yoh,w — 2 x Yh 2 o provides the best fit to the observed N(OH)/N(H 2 0) ratios. The same 
model also provides predicted 12 CO and atomic C columns and line intensities of the 12 CO 
J=l-0, J=3-2 and [CI] 609 fim lines, all of which were observed in the W49 components. 
The models overpredict the columns and the intensities of the CO J=l-0 and [CI] lines by a 
factor ~ 3 which may be explained by dilution of the emission lines in the beam. The J=3-2 
line, however, is predicted to be a factor ~ 2 smaller than is observed; a better fit for this 
line is found at higher densities, but makes the fit of all the other observed lines poorer. We 
note that Bensch et al (2003) had similar problems in fitting a PDR model to the diffuse 
cloud CO and [CI] results. 

Comparing our H 2 column to the model Ayt, we find a column-avearged H 2 abun- 
dance of 3 x 10~ 9 with respect to H nuclei. The observational papers cited above quote 
H 2 0/H 2 abundance ratios > 10~ 7 . These high abundances, however, were derived from the 
observed CO columns by multiplying the CO columns by 10 4 to obtain H 2 columns. In our 
models, we find H 2 columns of ~ 4 x 10 21 cm -2 and H 2 0/H 2 abundance ratios of ~ 10~ 8 . 
The H 2 /CO column-averaged abundance ratio in the model is 1.3 x 10 5 . In these translucent 
clouds, CO is photodissociated much more than H2, which self-shields much more effec- 
tively. Considerable H 2 exists in regions where the gas phase carbons is mostly C + . These 
translucent clouds are examples of "dark H 2 " which is not traced by CO. 
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4.2. Dense, Opaque Clouds 

4-2.1. Ay Thresholds for Ice Formation 

Whittet et al (2001) have presented plots of the water ice columns versus cloud extinction 
in which a clear threshold Ay for water ice is inferred. In our model we identify that threshold 
Ay with Ayf, since the ice very rapidly builds up with increasing Ay once the grain has 
a monolayer of ice. Whittet et al (2001), Tanaka et al (1990), and Williams, Hartquist & 
Whittet (1992) find that Ayf ~ 3 for Taurus and ~ 6 for p Oph. Assuming the radiation 
field in Taurus is Go ~ 1 an d that the gas density is 10 3 cm -3 , we predict from our analytic 
model with Y = 3 x 10~ 3 (i.e., Y H . 2 o = 10~ 3 , recall that Y is the tota l yield of both H 2 Q 



and of OH+H being ejected from water ice) that Ayf = 3. For p Oph, iLiseau et al.l (119991 ) 
and Kamegai et al (2003) have estimated G ~ 100 — 1000 as an average value in the SWAS 
beam. Assuming n = 10 3 — 10 4 cm -3 in the ambient gas averaged over the beam, we predict 
Ayf ~ 6. Thus, our model does a good job of matching these thresholds. Note that the 
threshold Ayf in our model mainly depends on hx(GoY/n) and so to the extent that we know 
Go/n, we can constrain Y. We therefore see that Y ~ 3 x 10 -3 is a reasonable choice for 
interstellar ice grains I 1 " 3 ! 



4.2.2. B68 

B68 is one of the best-studied, nearby, isolated molecular clouds (Alves et al. 2001, 
Zucconi et al. 2001). It lies at ~ 160 pc, with a mean density of n ~ 10 5 cm -3 , an A 
through the cloud center of ~ 30, and a central temperature of 7 K. Extinction and dust 
continuum maps have recently provided a detailed view of the density and temperature 
structure of the cloud (Alves et al. 2001, Zucconi et al. 2001). Previo us models of the CO 



v 



lines and the exterior temperature of the cloud suggest G ~ 0.25 — 1 (IBergin et al. II2002I ). 
The angular size of the cloud is ~ 1.3', smaller than the SWAS beamsize. 

SWAS has performed long integrations on B68 searching for the 557 GHz H2O line 
and the 487 GHz O2 line, setting stringent 3a upper limits on the line strengths (assuming 
thermal linewidths) from these undetected transitions of ^557 < 4 x 10~ 9 erg cm -2 s" 1 sr _1 
and ^487 ^, 6 x 10~ 9 erg cm" 2 s _1 sr _1 . Utilizing the known distribution of density n with 
Ay, and assuming Go = 0.5 and our standard values of Y and o~h, we have modeled the H2O 
and O2 abundance profiles and the expected observed line strengths in B68. To compare 



See also Nguyen et al (2002) for another model of the ice distribution in Taurus 
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with the observed intensities, we calculate the intensity expected in the model and then 
dilute the intensity by the cloud area divided by the SWAS beam area. Averaged over 
the SWAS beam, the predicted line intensities of ^557 ~ 4 x 10" 11 erg cm" 2 s" 1 sr" 1 and 
J 487 ~ 1 x 10" 12 erg cm" 2 s" 1 sr" 1 lie well below the upper limits set by SWAS. We predict 
a non-beam-diluted line intensity J1113 ~ 10" 8 erg cm" 2 s" 1 sr" 1 , which is not detectable by 
Herschel. Because the 557 GHz and 487 GHz line intensities are sensitive to Go (via the gas 
and dust temperatures), regions like B68 with low FUV fields are not promising targets for 
H 2 and 2 detections. 

Another interesting effect in our modeling of B68 arises because our models incorporate 
the observed increase of density with depth Ay. From Figure 6 one sees that this reduces 
the width A Ay of the H 2 plateau since the lower density surface gas means that A V f is 
larger, whereas the higher density gas in the interior makes Ayd smaller. However, this effect 
is not large, since our B68 model has n = 1.3 x 10 5 cm -3 at Ay = 1, n = 1.7 x 10 5 cm -3 
at Ay = 2, and n = 2.4 x 10 5 cm" 3 at Ay = 3 (roughly the extent of the H 2 plateau). 
We find that our predicted column of H 2 in our density-varying model differs by only 10% 
from a constant density, n = 1.7 x 10 5 cm" 3 , model. 



4.2.3. NGC2024 

NGC2024 is a dense, n ~ 10 5 cm" 3 s tar forming molecul ar cloud in Orion illuminated by 



a relatively intense, Go ~ 10 4 , FUV field (iGianini et al.ll2000l ). Contrasted to B68, it provides 



a test of our models for elevated radiation fields (i.e., G > 500). SWAS has observed both 
the H 2 557 GHz line and the 2 487 GHz line in this source. We have reanalyzed the 
SWAS data from the archives and obtain a beam-averaged intensity of 3 x 10~ 7 erg cm" 2 
g -i gr -i f or 557 GHz line and a 3 sigma upper limit on the 487 GHz line of 4 x 10~ 8 erg 
cm" 2 s" 1 sr" 1 [see also Snell et al (2000) for the 557 GHz line and Goldsmith et al (2000) 
for the 487 GHz upper limit]. 

Our models (see Figs. 12 and 13) reproduce the H 2 557 GHz line with beam averaged 
gas densities n ~ 10 4,8 cm" 3 and Go ~ 10 4 , assuming that the ortho/para H 2 ratio is in 
LTE with the grain temperature and the ortho/para ratio of H 2 is the steady state value from 
the formulation by Burton et al (1992). Here, in our peak water plateau, the gas temperature 
is ~27 K, and the ortho/para ratio of H 2 is ~ 1/6. The grain temperature is ~ 42 K and the 
ortho/para ratio of H 2 is 2.7. In the model the 2 487 GHz line intensity is predicted to be 
4 x 10~ 9 erg cm" 2 s" 1 sr" 1 , well below SWAS detectability. The peak (plateau) abundances 
for H 2 and 2 in the model are (see Figs. 5 and 6) 10~ 7 and 3 x 10" 6 respectively. The 
column averaged abundances, assuming N ~ 4 x 10 22 cm" 2 or Ay ~ 20 through the cloud 



-40- 



( IGianini et al.ll2000l . based on CO column densities greater than fewxlO cm over a ~ 4' 
region) is x{H 2 0) = N(H 2 0)/N ~ 1.1 x 10" 8 , and x(0 2 ) ~ 1.7 x 10" 7 . We conclude that 
the models also provide good fits to NGC2024, using reasonable values of n and Go, and 
predict peak abundances that are approximately 10 times greater than the column averaged 
abundances. We predict Im 3 ~ 3 x 10~ 7 erg cm -2 s" 1 sr -1 , which is detectable by the HIFI 
instrument on the Herschel Observatory. 



4-2.4- 2 Observed by Odin in p Oph 

Larsson et al (2007) have reported the detection of the ground state transition of 2 at 
119 GHz with the space submillimeter telescope Odin. Note that this is a different transition 
than the O2 transition observed by SWAS, and is at a much longer wavelength. Therefore, 
although the Odin telescope has twice the diameter of SWAS, the Odin beam size at 119 
GHz (~ 10') is more than twice that of the SWAS beam at 487 GHz. At the distance of 
p Oph, the Odin beam corresponds to 0.4 pc. The beam was centered on a 450 pm peak, 
and the submillimeter dust emission and the CO isotopic rotational emission suggest an H2 
column of 2 x 10 22 cm -2 (Ay ~ 20). A temperature of 30 K has been obtained for both 
the gas (e.g., Loren et al 1990) and the dust (Ristorcelli et al, in preparation) in this region. 
Using these parameters, and assuming the O2 to be optically thin and in LTE, Larsson et al 
find a column N(0 2 ) ~ 1 x 10 15 cm -2 for an abundance relative to H 2 of ~ 5 x 10~ 8 but an 
abundance in our units (relative to H nuclei) of x(0 2 ) = 2.5 x 10 -8 . The authors acknowledge 
a possible uncertainty of at least a factor of 2 in their estimate of the O2 abundance. 

The p Oph region observed is illuminated on the back side by a B2 V star, with the cool 
dense cloud core sitting in front of the photodissociation region. The relatively high dust 
temperature, and the likely di stance of the B star f r om the cloud, sugges t Go ~ 100 — 1000 



shining on the cloud surface (jCesarsky et al.l 1 19761 ; iKamegai et al.ll2003l ). The 0.4 pc size 



coupled with the observed total H 2 column suggest n ~ 10 5 cm -3 in this dense region. 
Setting Go = 300 and n = 10 5 cm -3 , and using standard parameters, we obtain a column 
N(0 2 ) = 1.4 x 10 15 cm" 2 . In our model, the gas and dust temperatures at the 2 abundance 
peak are ~ 20 K. We predict (see Figs. 12 and 13) an H 2 557 GHz intensity of J557 = 10" 7 
erg cm -2 s" 1 sr" 1 , and a 2 487 GHz intensity of I^-j = 2 x 10" 9 erg cm -2 s" 1 sr" 1 from 
this region. We have analyzed SWAS archival data on this source and find a measured H 2 
557 GHz intensity of 1.3 x 10~ 7 erg cm" 2 s" 1 sr" 1 and a 3a upper limit of 2 x 10" 8 erg cm -2 
s" 1 sr" 1 for the 2 487 GHz line. We predict /1113 ~ 3 x 10" 7 erg cm -2 s" 1 sr" 1 , which 
may be marginally detectable by the HIFI instrument on Herschel. In summary, our model 
successfully fits the H 2 and 2 observations of p Oph, using reasonable parameters for this 
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region. 

4-2.5. Upper Limits on 2 Observed by Odin 

Pagani et al (2003) report 3a upper limits on N{02) in a number of sources. As examples 
of sources with low incident G , the upper limits for TMC1 and L183 (L134N) are 7x 10 14 and 
1.1 x 10 15 cm -2 , respectively. Although the Odin beam is large, these sources are extended 
enough that beam dilution should not be a significant factor. In addition, although much 
of the Taurus cloud has relatively low total Ay, these regions appear to have sufficiently 
high Ay > 10 to encompass the 2 plateau. Using an estimated G ~ 1 and n ~ 10 4 cm -3 , 
our models predict a column of N(0 2 ) ~ 4 x 10 14 cm" 2 (see Figure 10) on both sides of 
the cloud, for a total column of ~ 8 x 10 14 cm" 2 . This prediction lies very close to the 
observed upper limits. If future observations or analysis drive this upper limit down below 
our predicted value, then either beam dilution or a lower value of Yh 2 o need to be invoked 
to reconcile observations with the models. Note that the upper limits on N(0 2 ) set upper 
limits on Y H , 2 o < 10" 3 . 

More challenging to our model are the upper limits quoted for sources with relatively 
high Go- For example, Pagani et al report 3a upper limits of N(0 2 ) < 1.9 x 10 15 cm" 2 for 
Orion A. This column is averaged over the 9' beam of Odin, which corresponds to a linear 
extent (diameter) of about 1.4 pc at the distance of Orion. If there were no extinction from 
the source of the FUV (primarily Q 1 C in the Trapezium), the FUV flux at a distance of 0.7 
pc from the star corresponds to Go ~ 4 x 10 3 . Taking this as the average Go incident on 
the molecular cloud in the beam, our model predicts N(0 2 ) ~ 10 16 cm" 2 (see Figure 10), 
significantly higher than the observed upper limit. We offer two possible explanations for this 
discrepancy. First, Q 1 C may be inside a cavity in the cloud carved out by the champagne flow 
created by the Trapezium. In this case, extinction by the intervening cloud would diminish 
G in the outer regions of the beam, and thereby diminish the beam average column of 2 . 
A second possibility arises from the fact that the enhancement of the 2 column at high G , 
as discussed in §3.2, is caused by the thermal evaporation of O atoms from grain surfaces 
before they can react to form water ice on the surface. We have adopted a binding energy of 
O atoms to the grain surface of Eo/k = 800 K (see Table 1). However, this binding energy 
is uncertain and in our view may be low, given that O is a radical. A modest increase in this 
binding energy would result in a much higher critical Go where the enhancement of N(0 2 ) 
occurs, since the grain temperature is a weak function of Go (T gr oc Gq 2 ) and thus the critical 
G goes as Eq. In other words, increasing E /k to 1200 K would increase the critical G 
from about 500 to about 4000, and provide a solution consistent with the observations. 
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Sandqvist et al (2008) use Odin observations to place an upper limit of N(0 2 ) < 6 x 10 16 
cm -2 toward SgrA. As noted in that paper, this is consistent with our model clouds for any 
incident G and any density n for standard local gas phase elemental abundances of O and 
C. 



5. SUMMARY, DISCUSSION AND CONCLUSIONS 

We have constructed both a detailed numerical code and simple analytical equations to 
follow the abundances of the main oxygen-bearing species as a function of depth in molecular 
clouds. We take clouds of constant density n, illuminated by an external FUV (6 eV< 
hv < 13.6 eV) field characterized by the unitless parameter G normalized to the local 
(Milky Way) interstellar FUV field. We explore the range 10 3 cm~ 3 < n < 10 6 cm -3 and 
1 < Go < 10 5 . Our models incorporate thermal balance, gas phase chemistry, a simple 
grain surface chemistry that includes the formation of H 2 , H 2 0, and CH 4 on grain surfaces, 
the sticking of gas phase species to grain surfaces, and various desorption processes such as 
photodesorption, thermal desorption, and cosmic ray desorption. 

The abundances of chemical species vary considerably as a function of depth into a 
molecular cloud. At the surface, molecules are photodissociated and the gas is primarily 
atomic (e.g., O) or, for species with ionization potentials less than 13.6 eV, singly ionized 
(e.g., C + ). Here, although dust grains may be cold enough to prevent rapid thermal desorp- 
tion, ices do not form because photodesorption keeps the refractory dust clear of ice mantles 
(the same is true of dust in the diffuse ISM). Very deep in the cloud, again assuming grains 
cold enough to prevent thermal desorption, photodesorption is negligible and ices form, 
causing gas phase abundances of condensibles to plummet. Here, time dependent effects 
(generally associated with cosmic rays) dominate the gas phase chemistry due to the time 
required to deplete condensibles onto grains, to desorb CO ice by cosmic rays, and to de- 
stroy gas-phase CO (by He + produced by cosmic rays). These processes are slow, and have 
timescales comparable to, or longer than, typical molecular cloud lifetimes. At intermediate 
depths, photodesorption of ices supply gas phase elements, and the partial shielding of the 
photodissociating external FUV flux allows gas phase molecular abundances to peak. 

We have focused in this paper on oxygen chemistry, and followed in particular the gas 
phase species O, OH, 2 , and H 2 as well as water ice. Our answer to the longstanding 
astrochemical question, 'Where is elemental O in molecular clouds?", is that oxygen not 
in CO is primarily gas phase atomic O to a depth Ay/ (see Eq. 23), and is water ice at 
larger depths. The freezeout depth Ay/ where significant water ice forms is proportional to 
ln(G /n), and is typically 3-6. Refractory silicate grains contain a modest abundance of O 
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(~ 10~ 4 ), but the bulk is likely gas phase O and CO along with CO ice and H2O ice. For 
example, in our standard model, integrating to Ay = 10, the O not in refractory grains is 
divided 56% H 2 0-ice, 25% gas phase O, 10% CO-ice, and 8% gas phase CO. 

Our models also give an indication of the origin of water in molecular clouds. In the 
steady state standard model, at Ayf the formation rate of gas phase H2O is 2% via gas 
phase reactions and 98% by photodesorption of H 2 that has been formed on grains. At 
A Vd the percentages are 30% and 70%, respectively; midway in the plateau, they are 8% 
and 92%, respectively! 1 " 4 ! These percentages depend on G and n, but as long as G < 500, 
these numbers are representative. At higher Go > 500, where O atoms thermally desorb 
from grains before forming OH and eventually water ice, most of the water production is 
via gas phase reactions. The production of water ice is about 3 times greater than the 
production of gas phase water by photodesorption of water ice, because 2/3rds of the water 
ice is photodesorbed as OH. Thus, one needs to modify these numbers if one wants to 
estimate the formation of water in any form in a cloud: for example, in the standard case 
midway in the plateau the percentages are 3% by gas phase reactions and 97% by grain 
surface reactions. 

We plan a follow-up paper to examine carbon chemistry. We note that Figure 14 shows 
that for all densities considered and by t > 2 x 10 6 years, gas phase CO is largely absent 
in cloud interiors. This is in rough agreement with observations that show CO depletions 
in dense cloud interiors. Thus, CO is not always a volume tracer. However, the sequel 
paper will include more grain surface carbon chemistry such as the formation of CO2 and 
methanol on grain surfaces. In addition, we will examine the effects of initial conditions 
on time-dependent chemistry, and compare observations with a wide variety of models with 
varying cosmic ray desorption rates and gas density profiles through the cloud. 

For Go ^,500, the abundances of gas phase OH, O2 and H 2 peak and plateau at 
values ~ 10~ 7 between A v f and Ay a (see, Eqs. 26, 31, and 33). These abundances are 
independent of n and Go because both their formation rate (ultimately from photodesorption 
of H 2 ice) and their destruction rate (from photodissociation) both depend on the product 
nGo- However, since their formation is linked with H 2 ice photodesorption, these plateau 
abundances do increase with increasing ice photodesorption yield, Y H2 o and Yoh,w, an d also 
with the grain cross sectional area per H nucleus an (see same equations). 

For higher values of Go > 500, the chemistry undergoes a change, as O atoms no longer 



1 A warning to future modelers: these percentages are not trivial to compute since there is a great deal 
of cycling of H2O and OH to H30 + (via reaction with Hj) and then back to H2O and OH (via dissociative 
electronic recombination). This cycling gives, of course, no net production of these molecules. 
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form H 2 ice on grain surfaces. Here, the grains are warm enough ( > 20 K) that O atoms 
thermally desorb before reacting with H atoms. This pushes the freezeout to greater depths 
and enhances the abundances of gas phase O, OH, 2 and H 2 by keeping more of the 
elemental oxygen in the gas phase. Because 2 formation goes as the product of O and OH 
abundances, 2 is the most sensitive to this change and its local abundance reaches 10~ 5 , 
while its column climbs to ~ 2 x 10 16 cm" 2 . We note that the observations of Pagani et al 
(2003) may suggest a higher binding energy for O on grains than adopted, and in this case 
the critical G where this occurs may be significantly higher than 500. 

Deep in the cloud interior at Ay » Ayf, the gas phase abundances drop as ices become 
dominant, the exact values of the abundances depending on cloud lifetimes and the cosmic 
ray desorption rates assumed (see Figures 14 and 15). Here, time dependent chemistry 
dominates. 

Poelman, Spaans, & Tielens (2007) have suggested that the SWAS results do not rule 
out substantial gas phase H 2 abundances and columns in the centers of clouds, because 
of optical depth effects which "hide" the H 2 in the centers. However, as this saturation 
effect occurs, the line center antenna temperature approaches the gas temperature of the 
emitting region. In our model, this gas temperature is typically 10-20 K. Even in the cold 
interior of clouds, one would expect temperatures of 5-10 K. The SWAS observed antenna 
temperatures in regions of detected H 2 are of order Ta < 1 K. This result alone suggests 
the H 2 is in the "effectively optically thin" regime. Our model does include an escape 
probability formalism for the H 2 lines. We find, in agreement with Poelman et al, that 
significant collisional de-excitation and saturation of the line occurs when Tn/n cr > 1, where 
r is the line center optical depth and n cr is the critical density (~ 10 s cm" 3 for the 557 GHz 
line). Thus, the requirement for 557 GHz line saturation is roughly 



N(H 2 Q) 



n 

n, 



> 1. 



(42) 



cr 



10 15 cm 2 and n < 10 6 cm 3 , we are almost always within 

10 6 cm" 3 . 



_ 6 x 10 12 cm" 2 

Since our models have N(H 2 0) ~ 

the effectively optically thin regime, except perhaps for the highest densities n 
Furthermore, if the H 2 columns and densities were high enough that line saturation oc- 
curred, the line fluxes would exceed those observed by SWAS. We argue that SWAS obser- 
vations indicate effectively thin 557 GHz emission, and that there is no evidence for hidden 
H 2 in the centers of clouds. 

Bergin et al (2000) listed ingredients for a successful astrochemical model of a molecular 
cloud. Slightly updated, they included: (1) low column-averaged abundances ( < 10" 7 ) of 
gas phase 2 , (2) low column-averaged abundances (<3 x 10" 8 ) of gas phase H 2 0, (3) an 
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explanation of the apparently higher abundances of H 2 when it is observed in absorption in 
translucent clouds, (4) column-averaged water ice abundances of ~ 10~ 4 , and (5) existence 
of complex carbon-bearing species in the gas phase in some regions. As noted above, the 
model presented in this paper fulfills all of these criteria. 

In addition, we have successfully applied, in a simple way, these models to observations 
of the Ay threshold for the onset of water ice, the strict upper limits for the gas phase H 2 
abundance in B68 made by SWAS, the possible detection of 2 in p Oph made by Odin, the 
detection of H 2 and upper limits for 2 made by SWAS for NGC 2024, and upper limits 
of the 2 abundance set by Odin. 

These successes imply that the proper modeling of the interstellar chemistry of molecular 
clouds must include as ingredients: photodissociation, freezeout on grains, grain surface 
chemistry, desorption processes, and, in the very opaque central regions, time-dependent 
chemistry. At the surface photodissociation and photodesorption dominate. Deep in the 
cloud, dust vacuums the condensibles and slow cosmic ray processes are important, requiring 
time-dependent chemistry. At intermediate depths, gas phase molecular species such as 2 
and H 2 peak in abundance, and provide the source of their submillimeter emission. 

We would like to thank M. Baragiola, U. Gorti, R. Liseau, L. Pagani, and E. van 
Dishoeck for many useful discussions. We also thank the referee for a thorough and construc- 
tive report. We gratefully acknowledge the financial support of NASA grant NNG06GB30G 
from the Long Term Space Astrophysics (LTSA) Research Program. 
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Table 1. Adsorption Energies, Photodesorption Yields, and Cosmic- Ray Desorption Rates 

of Surface Monolayer Atoms and Molecules 



Species Adsorption Energy 3, Photodesorption Yield Cosmic-Ray Desorption Rate 





[K] 


[UV Photon]^ 1 


[mol 1 s *] 





800 


10" 4 




OH 


1300 


10~ 3 




H 2 


4800 b 


10~ 3 as H 2 


4.4 x 10- 17 c 






2 x 10~ 3 as OH 




o 2 


1200 


10^3 


9.1 x 10- 15 d 


c 


800 


10-3 




CH 


650 


10-3 




CH 2 


960 


10-3 




CH 3 


1200 


10-3 




CH 4 


1100 b 


10-3 


2.7 x 10- 15 d 


CO 


960 b 


10-3 


6.0 x 10- 13 c 


s 


1100 


10-3 




Si 


2700 


10-3 




Fc 


4200 


10"3 





All adsorption energies are from lHasegawa fc Herbstl (119931 ) unless otherwise noted. 



Aikawa et all (119961 ) 



Herbst fc Cuppenl (120061 ) . Bringa & Johnson (2004) 



Hasegawa fc Herbstl (119931 ) 
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Table 2. Standard Model Parameters 



Parameter 


Symbol 


Value 


H nucleus density 


n 


10 4 cm" 3 


FUV field 


G 


100 


Oxygen abundance 11 


x 


3.2 x 10" 4 


Carbon abundance 3. 


x c 


1.4 x 10" 4 


Sulfur abundance 3. 


x s 


2.8 x 10~ 5 


Silicon abundance 3 


X S i 


1.7 x 10~ 6 


Magnesium abundance 3 


%Mg 


1.1 x 10~ 6 


Iron abundance 3 


X Fe 


1.7 x 10~ 7 


Total grain cross section per H nucleus 




2 x 10~ 21 cm 2 



3 Abundances of nuclei in the gas phase at the cloud surface, relative to H nuclei; Savage 
& Sembach (1996) 
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Table 3. Symbols Used in the Text 



Symbol Definition 



a grain radius 

Ay visual extinction from surface of the cloud 

A V f visual extinction where there is onset of water ice freezeout, monolayer of ice forms 

A Vd visual extinction where the gas phase water abundance drops 

AA V width of the plateau of gas phase H 2 or 2 , or (A V d — A V f) 

Ei binding energy of species % to grain surface 

f Q , (/ p ) fraction of H 2 in the ortho(para) state 

f S) i fraction of surface sites occupied by species % 

Fpuv fl ux of FUV photons at the cloud surface 

F Q 10 8 photons cm -2 s -1 , the flux of FUV photons corresponding to G — 1 

F pd i flux of molecules of species i photodesorbing from ice surface 

F t d,i flux of molecules of species i thermally desorbing from ice surface 

Go unitless parameter measuring incident FUV flux normalized to 1.6 x 10 -3 erg cm~ 2 s 

or to 10 8 photons cm -2 s -1 , roughly the local interstellar value 

7n e + net rate coefficient for destruction of H 2 by He + 

7 H + net rate coefficient for destruction of H 2 by Hjj" 

7o 2 rate coefficient for gas phase reaction of OH with O to form 2 

J 4 87 intensity of the 2 487 GHz transition, in erg cm -2 s _1 sr^ 1 

-Z557 intensity of the H 2 557 GHz transition, in erg cm~ 2 s _1 sr _1 

/1113 intensity of the H 2 1113 GHz transition, in erg cm~ 2 s _1 sr _1 

n gas phase hydrogen nucleus number density [~ n(H) + 2n(H 2 + n(H + )] 

n gr number density of grains 

V; L vibrational frequency of species i bound to a grain surface 

gas-phase column density of hydrogen nuclei 

Ni gas-phase column density of species % 

Nf column density of hydrogen nuclei for onset of water ice freezeout 

N s ^ number of adsorption sites per cm 2 on a grain surface (~ 10 15 cm -2 ) 

R t d,i thermal desorption rate of per atom or molecule of species i from a surface 

Rcr,% cosmic ray desorption rate per surface atom or molecule, i 

Rh 2 o 5.9 x 1CT 10 s _1 , the unshielded photodissociation rate per molecule of H 2 for G =l 

Ro 2 6.9 x 1CT 10 s -1 , the unshielded photodissociation rate per molecule of 2 for G =l 
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Table 3 — Continued 



Symbol Definition 



i?OH 3.5 x 1CT 10 s _1 , the unshielded photodissociation rate per molecule of OH for G =l 

oh grain cross sectional area per H nucleus 

<j gr cross sectional area of a single grain 

t cr>i timescale for cosmic ray desorption of adsorbed ice species % with abundance x(i) 
t e vap,u 2 o time for an adsorbed H 2 molecule to thermally evaporate from a grain surface 

t e vap,o time for an adsorbed atom to thermally evaporate from a grain surface 

tfj timescale for species % to hit (and stick to) a grain 

ty timescale for a grain of size a > 100 A to absorb an FUV photon at the cloud surface 

tjj timescale for a grain of size a > 100 A to absorb an FUV photon at A V j 

try-des timescale to photodesorb an O atom from a grain with covering factor / S; o 

tgrj time for a given grain to be struck by a gas phase atom or molecule of species % 

t r timescale for a grain to radiate much of its thermal energy 

T ave the average of the gas T at Ay/ and at the Ay where the H 2 abundance peaks 

T m ax peak temperature of a grain after absorbing one FUV photon 

^300 T/300 K, where T is the gas temperature 

x(i) abundance of species % relative to H nuclei 

x Q abundance of ortho-H 2 

x p abundance of para-H 2 

x p i t i abundance of species i in the plateau region between Ay/ and Ayd 

Yi photodesorption yield of species i from a grain surface 

V Vh 2 o + V h>, the total yield of H 2 and OH being desorbed from water ice surface 

Z(T) partition function for 2 
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Fig. 1. — Standard ion- neutral chemistry leading to the formation of H 2 and 2 in shielded 
regions of the ISM. The chemistry is often initiated by cosmic ray ionization of H2. The 
resultant reacts with H2 to form hLt. An alternate route starts with the cosmic ray 
ionization of hydrogen and the resultant change exchange of H + with O to form + . 



OH H 2 



Fig. 2. — Grain-surface oxygen chemical network included in this modeling effort. Note that 
in addition we include grain surface carbon chemistry (see text) and allow the sticking of all 
species onto grain surfaces and various desorption mechanisms. The grain surface oxygen 
reaction network begins with O sticking to grains, followed by accretion of H to form OH ice 
and then H 2 O ice . Here, "d" stands for the three desorption processes: thermal, photo-, 
and cosmic ray. Near the cloud surface, FUV photons photodesorb ices, leading to some 
gas-phase H2O and O2. In addition, for high FUV fields which warm the grains, thermal 
desorption of the weakly bound O atoms can be important both near the surface and deep 
in the cloud. Finally, deep in the cloud, cosmic ray desorption is important. 



2 4 6 8 10 

A v 



Fig. 3. — Oxygen chemistry for the standard (n = 10 4 cm -3 , Go = 100) steady-state PDR 
model, including freeze- out, grain-surface reactions, thermal desorption, photodesorption, 
and cosmic ray desorption. Model parameters are those given in Table 1. The abundances of 
major O-bearing species are shown as a function of visual extinction from the cloud surface, 

Ay. 



2 4 6 8 10 

A v 



Fig. 4. — Carbon chemistry for the standard (n = 10 4 cm -3 , Go = 100) steady-state PDR 
model, including freeze- out, grain-surface reactions, thermal desorption, photodesorption 
and cosmic ray desorption. Model parameters are those given in Table 1. The abundances 
of major C-b earing species are shown as a function of visual extinction from the cloud surface, 
Ay. For Ay > 5, the gas phase CO abundance will not be as low as is shown here for the 
steady-state solution; time- dependent effects will keep the CO abundance elevated (see §3.7). 
In addition, the limited chemistry in the steady state PDR code drives the high abundance 
of CH4 ice at high Ay; the time dependent code, which has more extensive chemistry, finds 
a mixture of CO and CH 4 ice for the carbon-bearing ices in this case. 
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Fig. 5. — H2O, O2 and H20j ce abundances for a cloud with n = 10 4 cm -3 but with a variety 
of FUV field strengths incident on the cloud surface. Results are shown for FUV fields Go=l, 
10 2 and 10 3 times the average interstellar field. H 2 Oj ce is dot-dash, H 2 is solid, and 2 is 
dotted lines. Higher Go drives curves to the right. Although the depth at which freezeout 
occurs is affected by G , the total H 2 column is not. The increase in the peak abundance 
of 2 seen for G = 1000 is caused by thermal desorption of atomic O from the warm grains, 
which suppresses H20j ce formation and keeps more elemental O in the gas phase. 



-59- 




Fig. 6. — Effect of changing the gas density. Results are shown for n = 10 3 , 10 4 and 10 5 cm -3 
and for the standard FUV field Go — 10 2 . The threshold Ay for water ice formation and 
the Ay where H 2 and O2 peak increase for increasing Gq/u. The peak abundances do not 
change with n. The labels on the arrrows refer to the log of the density n. 



-60- 




Fig. 7. — The gas temperature T as a function of Ay for a variety of n and Go with all 
other parameters standard. Dashed lines refer to gas density n = 10 4 cm -3 , and in order 
of increasing T at low Ay correspond to Go — 1, 10, 100, 10 3 , and 10 4 . Dotted lines refer 
to gas density n = 10 5 cm -3 , and in order of increasing T correspond to Go = 10 4 and 10 5 
(Orion PDR conditions). The solid portion of each curve refers to the region of the "water 
plateau" for G < 500, or to regions where the water abundance is at least 0.33 times its 
peak value for Go > 500. 
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Fig. 8. — Abundances of H 2 0, H 2 0i ce and 2 for the standard case (n = 10 4 cm -3 , G = 100) 
as a function of Ay and of the photodesorption yield of water ice to form gas-phase H 2 0, 
Yh 2 o Results are shown for the standard yield of 10~ 3 and for yields of 10~ 2 and 10~ 4 . 
H20i ce also photodesorbs OH with a rate two times 1h 2 o ( see text). The labels refer to the 
log of the yield. 
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Fig. 9. — Abundances of H 2 0, 2 and H 2 ice as a function of Ay and of the total grain 
cross section per H nucleus, er# . The higher curves for 2 and H 2 and the lower curves for 
H 2 0j ce refer to higher values of an- 
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Fig. 10. — The total columns of H 2 (left panel) and O2 (right panel) as a function of Go 
and n, assuming clouds with Ay greater than the plateau values, or roughly Ay > 3 — 6. The 
increase of the column of H 2 and O2 at Go ~ 500 is caused by the thermal desorption of O 
atoms from grain surfaces, which suppresses the formation of water ice and enhances the gas 
phase abundance of elemental O. Note that since the column is produced at intermediate 
Ay, the column average abundance of H 2 and 2 will fall as Ay 1 for high Ay. Results are 
shown for the standard photoyield Yu 2 o — 10~ 3 for densities of 10 3 , 10 4 , and 10 5 cm -3 , and 
for Yh 2 o — 10 -4 for a density of 10 4 cm~ 3 . 
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Fig. 11. — The average (see text) gas temperature in the water plateau region or the peak 
water abundance region for a range of densities and FUV field strengths incident on the 
cloud surface. Note that the temperature only increases by a factor of ~ 3 while the incident 
FUV flux G increases by 10 5 . The plateau retreats deeper into the cloud as Go increases, 
which suppresses a more dramatic increase in T. 



-65 - 




Fig. 12. — H 2 intensities for a range of densities and FUV field strengths incident on the 
cloud surface. Contours are labeled with the logarithm of intensity, measured in units of 
erg s _1 cm" 2 sr _1 . Left Panel: Intensity of the ortho-H 2 groud state transition li,o — lo,i at 
557 GHz. To convert to units of Kkms -1 , multiply contour levels by 5.6 x 10 6 ; Right Panel- 
Intensities of the para-H 2 ground state transition — ,o at 1113 GHz. To convert to 
units of Kkms -1 , multiply contour levels by 7.0 x 10 5 . 
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logi [ I (0 8 487 GHz)/erg s" 1 cm" 2 sr" 1 ] 




log 10 [n/cm- 3 ] 

Fig. 13. — The O2 487 GHz intensity for a variety of densities and FUV field strengths 
incident on the cloud surface. Contours are labeled with the logarithm of intensity, measured 
in units of erg s _1 cm -2 sr" 1 . To convert to units of Kkms -1 , multiply contour levels by 
8.4 x 10 6 . Note the relative constancy of the O2 intensity for Go < 100, since the column 
and temperature of the O2 plateau region stay quite constant. 
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Fig. 14. — Abundances of H2O (solid), CO (dot-dash) and O2 (dotted) at Ay = 10 at a 
cloud age of 2 x 10 6 years, for a range of CO cosmic-ray desorption rates. Our standard CO 
cosmic-ray rate is given by k(CO cr )o and is given in Table 1. Results are shown for rates 
from 1CT 2 to 10 2 times the standard rate, for Go = 100 and densities % = 10 3 , 10 4 , and 
10 5 cm -3 . The higher densities have lower abundances due to freezeout of gas species. Note 
that the H 2 and O2 abundances are only significant (i.e., comparable to the plateau value), 
for low densities. The timescale t cr< co for cosmic ray desorption of CO ice (see Eq. 8) is 
given at the top. 
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Fig. 15. — Same as Fig. 14 but for an age of of 2 x 10 7 years. Note that the H2O and 
O2 abundances are never significant compared to their plateau values. If the cosmic ray 
desorption timescale for CO ice, t cr> co is l° n g; the H 2 and 2 abundances are low because 
the frozen CO is not liberated to the gas to provide elemental gas phase O to form them. If 
the timescale is very short, the CO ice is desorbed quickly and the gas has time to deplete 
the elemental gas phase O as water ice on the grain surfaces. 
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Fig. 16. — Ratio of carbon nuclei to oxygen nuclei in the gas phase for clouds with a variety 
of densities and FUV field strengths. Results are shown for a cloud age of 2xl0 6 years. If 
the C/O ratio rises above unity, then interesting organic chemistry ensues (see text). Note 
that the C/O ratio only rises above unity for one case (the high density, n = 10 5 cm -3 case). 
However, in this case the time to deplete the gas phase carbon is very short so that the gas 
phase elemental C abundance is very low (bottom panel). Thus, the abundance of organic 
molecules will be very low. 
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Fig. 17. — Same as Figure 16 but for cloud age of 2xl0 7 years. Note that C/O rises above 
unity for all cases except the case with high G = 10 3 . Thus, in all cases but this one, 
complex carbon chemistry is initiated, since not all C will be tied up in CO. The bottom 
panel shows which clouds will have the highest abundance of these molecules (the cases with 
C/0> 1 and high elemental C abundance). The elemental C abundance rises with increasing 
G /n. 



